A linearized stabilizer formalism for systems of finite dimension 
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The stabilizer formalism is a scheme, generalizing well-known techniques developed by Gottes- 
man [l| in the case of qubits, to efficiently simulate a class of transformations (stabilizer circuits, 
which include the quantum Fourier transform and highly entangling operations) on standard basis 
states of d-dimensional qudits. To determine the state of a simulated system, existing treatments 
involve the computation of cumulative phase factors which involve quadratic dependencies. We 
present a simple formalism in which Pauli operators are represented using displacement operators 
in discrete phase space, expressing the evolution of the state via linear transformations modulo 
D ^ 2d. We thus obtain a simple proof that simulating stabilizer circuits on n qudits, involving any 
constant number of measurement rounds, is complete for the complexity class coModdL and may be 
simulated by 0(log(n) 2 )-depth boolean circuits for any constant d ^ 2. 
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I. INTRODUCTION 

Efficiently simulating the evolution of quantum states 
is a problem of central importance to modern physics. 
It is generally considered to be a difficult problem: this 
is formalized by the conjecture that the complexity class 
P, representing decision problems efficiently and deter- 
ministically solvable by conventional computers, does not 
contain BQP, the class of decision problems solvable effi- 
ciently and with bounded error by using local many-body 
interactions acting on a standard basis state. 

However, some transformations of quantum states are 
indeed easy to classically simulate. The stabilizer for- 
malism is one technique to do so, which may be applied 
for some evolutions of finite-dimensional systems. Devel- 
oped by Gottesman 1] to study quantum error correc- 
tion on qubits, the stabilizer formalism describes how to 
simulate stabilizer circuits: CP maps obtained by com- 
posing Clifford group operations (unitary maps which 
preserve the group generated by Pauli spin operators 1, 
a x , <j y , <j z ), Pauli observable measurements, and Clif- 
ford operations controlled directly on measurement out- 
comes 0- This serves as a foundation for other efficient 



simulation techniques 
quantum computing 



and for measurement-based 
It has since been extended 
in similar ways by Gottesman Q and by Hostens, De- 
haene, and de Moor [9j to systems of qudits of dimension 
d ^ 2. This generalization of the stabilizer formalism 
can be used to simulate circuits involving the quantum 
Fourier transform (an important element of Shor's inte- 
ger factoring algorithm |l0() and certain arithmetic oper- 
ations on standard basis states (which are closely related 
to other key operations involved in Shor's algorithm). 

The stabilizer formalism is notable in part for the min- 
imal computational power it requires. Aaronson and 
Gottesman [ll[ show that for qubits, simulating stabi- 
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lizer circuits on computational basis states is feasible for 
the complexity class ©L, the subclass of P consisting of 
those problems which are reducible (using logarithmic 
workspace) to solving systems of linear equations mod- 
ulo 2 [l2j . However, in existing algorithms [0, [Hj], the 
exact state of the system is determined by scalar factors 
which depend quadratically on the exponents of Pauli op- 
erators cr®<Jz ; similar remarks apply to the formalism of 
Ref. [9] for d ^ 2. One may ask whether these quadratic 
dependencies are necessary for the efficient simulation of 
stabilizer circuits, or if they may be removed to obtain a 
formalism involving only linear algebra. 

This article presents a simple modification to the for- 
malisms of Refs. @, [H|, to simulate stabilizer circuits by 
solving systems of linear equations modulo D G {d, 2d}. 
Using discrete Weyl operators as a matrix basis, rather 
than the related operators o~%o~ z , we eliminate quadratic 
dependencies from simulations to enable Clifford group 
operations to be simulated entirely linearly. Further- 
more, we demonstrate a subgroup of Clifford group which 
forms a group representation (acting on density opera- 
tors) of the symplectic group modulo D, simplifying the 
connection to symplectic operators described in Ref. j^j. 
This allows us to obtain direct proofs, for any constant 
d 2, that simulating a constant number of measure- 
ments in a stabilizer circuit is complete for the complexity 
class coModdL C P defined by Buntrock et al. [14| gen- 
eralizing ©L to arbitrary modulus d 2. Consequently, 
all such circuits can be simulated by 0(log(n) 2 )-depth 
boolean circuits. 



Beyond the complexity theoretic results of this article, 
the formalism that we present here is practically useful 
for simulations of stabilizer circuits. This (as well as fur- 
ther technical results on Pauli stabilizer groups in the 
Appendices) should be useful to the study of error cor- 
rection [l5l.ll6| and measurement-based computation [T7j 
on higher-dimensional systems. 
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II. PRELIMINARIES 

Throughout, d ^ 2 is a constant. Let T~Ld — C d with 
standard basis vectors \q) for q £ {0, 1, . . . , d— 1}. When 
little confusion may result, we identify g with an element 
of Zd (the ring of integers mod d) . 

Our formalism differs slightly from traditional descrip- 
tions of the stabilizer formalism, to achieve as uniform 
a description as possible while remaining consistent with 
conventions e.g. for the special case d = 2 for qubits. We 
first address a non-uniformity noted by Refs. [1, Q be- 
tween the cases of d even and d odd. In analogy to the 
case d = 2, we define unitary operators X and Z by 

X\q)=\q + 1) and Z \q) = e 2 ^' d \q) , (1) 

using addition mod d; then X and Z have order d. We 
also wish to have an operator Y which has order d such 
that XYZ = rl, where t 2 = e 2ll ^ d (generalizing the 
equation XYZ = il which holds for qubits). One may 
show that (X^Z^f = 1 for d odd, and (X^Z^) d = -1 for 
d even: the order of r must then depend non-trivially on 
the parity of d. We choose t = ^i^+^/d = (-^d^iv/d^ 
which is consistent with t — i for qubits [18| : then 
(tX^Z^ Y = 1 for all d. 

Definition I: r = c^^ 1 )^, let X and Z be defined 
as m Eq. ((T|), and Y = tX^ Z^ . The Pauli group over Z^ 
is the group V d = (X,Y,Z) = (t1,X,Z). The n-qudit 
Pauli group p® n is the group generated by n-fold tensor 
products of 1, A, Y, and Z. 

To reduce computations for phases in the stabilizer for- 
malism, we describe Pauli operators in terms of trans- 
lation operators in discrete phase space [l9|, [2(| , rather 
than directly as multiples of Z a X b : 

Definition II: A Weyl operator (on a single qudit) is 
an operator W a ^ = r~ ab Z a X b for some a, b £ Z . On n 
qudits, for vectors a, b S Z™, the Weyl operators are 



Wa,b = W aiM ® • • • <g> W 0nAl 

= T - ah (Z ai «) ■■■(g) Z Q ")(X bl 



XH (2) 



where a • b is the dot product. For block vectors v = 
a © b £ Z 2n , we write W v = Ifa.b- 

The Pauli operators {1,X,Y,Z} are all Weyl operators 
(for example, Y = W-i.^i); and the set of Weyl opera- 
tors is invariant under cyclic permutations of {X, Y, Z) in 
Definition |TI1 Weyl operators have convenient properties: 
define the symplectic inner product for v, w £ Z 2n by 



[v, w] = v T 02riW, where a 2ri 



On K 



(3) 



Note that [v, w] = — [w, v], and [v, v] =0. The following 
may be verified by considering actions on the standard 
basis, using the relation X b Z a = T -^z a X b for a,beZ: 



Lemma 1: For any v,w € Z 2n , we have fl^Ww — 
T K W ]VK V+W . The following properties then also hold: 

(i) W V W W = t 2 ^^W w W v . 

(ii) [W v , W w ] = if and only if [v, w] = (mod d). 
(Hi) Wl — Wt v for t e Z; in particular, W$ = W~ v . 
(iv) The order of W v divides d. 

These properties are a straightforward generalization of 
observations made by Appleby [T^] in his nearly iden- 
tical formalism for single qudit operators; they are also 
essentially what allow us to define the formalism of this 
article. They also suggest a convention of evaluating each 
vector modulo d: however, this leads to inconsistencies 
for d even. For instance, we have 



W ,i = Z°X 1 = r 2d Z d X 1 = -W d , 



(4) 



as Wd,i is non-zero, we then cannot equate Wo,i to Wd,i ■ 
Despite this, it is still possible to use modular arithmetic 
with the Weyl operators, as outlined in the next lemma: 

Lemma 2: For all v, w € Z 2 ™, we have W v oc Ww if and 
only if w = v + dx for some x e Z™ , in which case 



W w = (-l^+^^Wv. 



(5) 



In particular, W v — W w if v = w (mod d) for d odd, 
and if v = w (mod 2d) for d even. 

This follows from the action of WvWX on the standard 
basis. We may subsume the cases of d even and d odd into 
a single formalism by adopting the following convention: 



Convention: For d ^ 2, we define 



D 



the order of r 



d if d is odd, 
2d if d is even; 



(G) 



then the mapping v i-> W v is well-defined for v e Z 2 ^ 1 , 
and all associated arithmetic may be performed mod D. 

We frequently describe operations in terms of arithmetic 
modulo D — with a notable exception: as standard basis 
states \q) G Wd involve integers q £ {0, 1, . . . , d — 1} by 
definition, we require expressions involving q occurring as 
powers of r to be invariant under replacing q with q±d. 

Consider the Hilbert-Schmidt inner product (A, B) = 
Tr(A> B)/d n , renormalized to obtain (1,1) = 1. From 
the the effect of P £ , p® n on standard basis states, one 
may easily show Tr(P) ^ if and only if P oc 1. Then 
for P,Q£ Vf n , either P oc Q or (P, Q) = 0, so that: 

Lemma 3: The following are equivalent for v,w £ Z 2ra , 
for any system of n qudits of dimension d: 

(<) (W v ,Ww)^0; 

(ii) (Wv,W w )=±l; 

(Hi) Wv oc W w ; 

Furthermore, when these hold for d odd, W v = Ww 



(w) W V = ±W W ; 
(v) v = w (mod d). 
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From this, it is easy to show that a single Weyl operator 
Wp can only be generated by a collection of other Weyl 
operators W Vl , . . . , W ve if W p = ±W Vj for some operator 
W v . . By a simple dimension-counting argument, and us- 
ing the fact that two integers ^ a, b < d are equivalent 
mod d if and only if they are equal, we may show: 

Corollary 3a: The operators W v for v e {0, . . . , d- l} 2n 
form an orthonormal basis for the linear operators on 
Hf n with the inner product (A,B) = Tr(A^B)/d n . 

Note that the Weyl operators W v for v 6 Z 2 -^ are not 
independent for D even; this is an important technical 
point in the analysis to follow. 



III. SIMULATING UNITARY STABILIZER 
CIRCUITS 

We now describe a simple formalism, based on Weyl 
operators, for simulating an important subclass of stabi- 
lizer circuits: the class of unitary Clifford circuits, acting 
on standard basis states, and including at most one fi- 
nal measurement. We refer to these as unitary stabilizer 
circuits. 

Note that unitary stabilizer circuits exclude interme- 
diate measurements, or Clifford operations controlled by 
measurement outcomes. If we are only concerned with 
the distribution of a single measurement as output, we 
can still simulate circuits with classically controlled Pauli 
operations conditioned on measurement outcomes, by ap- 
plying the principle of deferred measurement (as we out- 
line in Section IIII C|) . It is also straightforward to ap- 
ply the techniques of this section to simulate classically 
controlled Clifford operations, provided the values of the 
controls do not arise from measurement outcomes. 

Unitary stabilizer circuits admit a simple representa- 
tion in terms of linear transformations modulo D, sim- 
plifying the individual operations used in the existing 
formalisms for d = 2 and d ^ 2 (in Refs. [l| and @). We 
use this simplified formalism to develop a more general 
stabilizer formalism which can describe the evolution of 
states under measurement, for arbitrary qudit dimension 
d, in Section ITVl 



A. Stabilizer tableaus for qudits 

We begin with a substitute for the so-called binary 
representation [ll|, EH of Pauli operators: 

Definition III: A Pauli vector on n qudits is a vector 
v e Z^ l+1 (which we represent as a column vector to 
enable transformations by left- multiplication) , which we 
decompose into blocks v = [(f) \ Vi v 2 ■ ■ ■ v 2n ] T ■ 

Every Pauli operator P with order ^ d is proportional 

to a Weyl operator by a power of r 2 . It suffices to show 
P d _ ( T h Wv y _ j if and only if T dh = X) and solve 



modulo D for h. Thus we may use Pauli vectors to rep- 
resent any Pauli operator of order at most d, using the 
correspondence 




As in the binary case, a stabilizer group is an abelian sub- 
group of Pauli operators which each have +l-eigenspaces. 
(Such operators P may be represented by Pauli vectors: 
we only have P d ^ 1 if d is even and P d = —1, in which 
case P has no +l-eigenspace.) Such a group has a joint 
+l-eigenspace: for a set § = {Si, . . . , Si} of generators, 
the projector onto the +l-eigenspace of S £ § is given by 

by expanding Tr(IIs 1 IIs 2 • • • Hs e ) as a sum of traces of 
Pauli operators, we find that the trace is non-zero. 

Definition IV: A stabilizer tableau Tg on n qudits is 
a matrix over Z^> with 2n + 1 rows, whose columns are 
Pauli vectors for the generators § = {Si, S%, . . . , St} of a 
stabilizer group. The first row </>g of the tableau is called 
the phase vector, and the rest Wg the Weyl block of the 
tableau. We may present Tg in terms of n x I blocks Xg 
and Zg , 



4>s 




<t>s 






Zg 















This definition corresponds closely to the notation in 
Refs. [HEldll f° r the binary case. The principal distinc- 
tion between this representation and the binary represen- 
tation is that the Zg and Xg blocks represent coefficients 
for Weyl operators, which differ from simple powers of Z 
and X operators on various qudits by scalar factors; and 
also that the coefficients are defined modulo D, rather 
than modulo d. 

The more important distinction is between Weyl oper- 
ators and simple powers of Z and X: while subtle, this 
difference is significant because of the product formula for 
Weyl operators in Lemma Q] Consider stabilizer tableaus 
satisfying a further constraint: 

Definition V: A stabilizer tableau is proper if the 
columns of its Weyl block Wg are orthogonal (mod D) 
with respect to the symplectic inner product; i.e. if 
Wjo-anWg ee (mod D). 

As the generators S € § of a stabilizer group commute, 
the columns of the Weyl block in a corresponding tableau 
are orthogonal mod d under the symplectic inner prod- 
uct by Lemma [1] proper tableaus simply satisfy the same 
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constraints mod D. Not all stabilizer groups can be rep- 
resented by proper tableaus (see Section IIV Al for a sim- 
ple counterexample); however, every standard basis state 
\QiQ2" ' On) S ^f"" nas a proper tableau of the form 



qi q_2 ■■■ q n 



0„ 



(10) 



We may show this by noting that \qf) £ 7id is the state 
stabilized by r~ 2qj Z acting on the j th qudit [2l|. Fur- 
thermore, every stabilizer group has a proper tableau for 
d = 2 or d odd; in Appendix IA 11 we show how to obtain 
such a tableau from an "improper" one. We may thus 
restrict our attention to proper tableaus in many cases. 

Given a proper stabilizer tableau, in which we repre- 
sent Pauli operators in terms of Weyl operators, we may 
represent recombinations of stabilizer generators using 
only linear operations on columns, with no further com- 
putation of phases: 



Lemma 4: For two Pauli vectors 4>\ © vi and <f>2 
drawn from a proper stabilizer tableau, we have 



-2,p 



-202 



W v 



-2(01+0 2 ) w 

y y V} 



+V 2 



v 2 



(11) 



This follows directly from Lemma [T] Compare this to 
multiplication of bare products of Z and X, in which 
(Z <8> Z){X ®X) = (ZX <g) ZX), but (X ® Z)(Z ®X) = 
— (ZX <E> ZX). In the usual binary representation used 
e.g. in Ref. [ll|, as well as the representation of Ref. [1] 
for arbitrary d ^ 2, supplemental phase-factors must be 
computed depending on which generators are being mul- 
tiplied; the phases depend quadratically on the exponents 
of X and Z involved in the operator. Our formalism ab- 
sorbs these quadratic dependencies into the calculus of 
Weyl operators, so that they do not appear elsewhere 
in computations involving products of generators. Thus, 
changes in the generating set of a stabilizer group can be 
represented using only linear transformations on a proper 
tableau. As we see in the next section, a representation 
using Weyl operators allows other transformations of the 
tableau to be effected linearly as well. 

For d even (e.g., for the case d = 2 corresponding to 
qubits), stabilizer tableaus are defined modulo D = 2d 
rather than mod d. This difference from Refs. [U,|Tl],[l3| is 
significant when using Weyl operators to represent Pauli 
operators, as illustrated in Eq. Q; and varying coeffi- 
cients of the tableau by ±d may not preserve the property 
of the tableau being proper. In Section ftV Al we describe 
a way in which the Weyl coefficients of a tableau may 
also be reduced mod d (though this involves techniques 
which accommodate improper tableaus). However, the 
phase vector of a stabilizer tableau may always be eval- 
uated modulo d, as differences of d in phase coefficients 
amount only to a scalar factor of r 2d = +1 for the oper- 
ator represented. 



B. The Clifford group on qudits 

Definition VI: The Clifford group on n qudits is the 
group of operators U £ U(Wf ™) such that UPW £ Vf n 
for any P £ Vf n ; we denote it by C n {d). 

1. Linear transformations of stabilizer tableaus 

For any U £ C n (d), the superoperator 3>(M) := UMW 
is an automorphism of matrix algebras; we may then re- 
duce the study of the Clifford group to its effect on a 
generating set {rl, Z x , Z 2 , . . . , Xi,X 2 , . . .} of V® n . Fur- 
thermore, $ maps rl to itself, and maps each Pauli op- 
erator of order d to another Pauli operator of order d. 
Thus we can characterize the effect of an arbitrary Clif- 
ford operator on n qudits by a (2n + 1) x (2n + 1) array 
over Zd (a "conjugation tableau"), 



1u 



h>2n 



Cu 



(12) 



whose columns are vectors representing the images of 
{t 2 1, We 1 , ■ ■ ■ , We 2 „} under $, in sequence. The effect 
of conjugating a set S of operators given by a tableau as 
in Eq. © by a Clifford unitary U can be evaluated by 
matrix multiplication modulo %u '■ 



(13) 



4>usu^ 




4>s 




= 7 V 












w s 



We demonstrate this by identifying a class of Clifford 
operators whose conjugation tableaus have coefficients 
hj = for each 1 ^ j ^ 2n: 

Definition VII: A symplectic Clifford operator U £ 
C n (d) is one for which UW^W is a Weyl operator for 
each 1 < j ^ 2n. We write the set of such operators as 

0C n (d). 

For U £ C n (d) with a tableau as in Eq. (|12l) . we may 
decompose U = U' W a2n h', by considering the effect of 
conjugating each of the Weyl operators Wg., by UW^lh , 
we may show that U' £ oC n (d) . We then prove Eq. (|T3]) 
by relating the Weyl blocks C\j< of U' £ aC n (d) to sym- 
plectic transformations modulo D: 

Definition VIII: A symplectic transformation of Z^ 1 is 
an operator C : Zj^ 1 —> such that C T G2 n C = <J2n , 
so that C preserves the symplectic inner product modulo 
D. We denote the group of such operators by Sp 2n (Z£>). 



Theorem 5: For any U £ oC n (d), there exists an oper- 
ator C £ Sp 2 „(Z£)) such that 



UW V U^ = W Cv , forallveZg 1 



(14) 



In particular, aC n {d) is a group. Similarly, for any C £ 
Sp 2n (Zr,), there is a U £ aC n {d) for which Eq. flU} holds. 
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Proof: Any C G Sp 2n (Z/j>) induces a superoperator $ 
defined on linear operators acting on Hf n , by extending 
$(Wv) = W Cv linearly [HI]. As C is invertible, so is 
$; then $ ® 1 is a matrix automorphism (where 1 is 
the identity superoperator on any auxiliary space). In 
particular, $ ® 1 maps projectors to projectors; thus it is 
positive, so that $ itself is completely positive. Finally, as 
$(1) = 1 and $(W V ) has trace zero for W v ^ 1, it follows 
that <& is trace-preserving. Then $ is an invertible CPTP 
map, so that $(M) = UMW for some J7 € oC n {d). 

Conversely, for U € oC(H® n ), the map M H- [/Af[/t 
preserves group commutators Wy^wW^W*, of Weyl op- 
erators within the Pauli group. This constrains the Weyl 
block Cu of Tc/ : for any 1 ^ h, j < 2n, if UW &h U ] = W Vh 
and UW^W = W Vj , we find that [v^,Vj] = [%,%] 
(mod d). Thus, the 2n x 2n array Cjj in fact represents 
a linear transformation of Z 2 ^ which is symplectic mod- 
ulo d. For d even, we then "lift" Cjj to obtain a matrix 
= Cu (mod which is symplectic mod 2d. (The 
procedure to do so is technical, and similar to the proce- 
dure for producing proper stabilizer tableaus for d prime 
or d odd: we defer this detail to Appendix IA2I ) □ 



The above Theorem essentially extends [19|, Lemma 2] to 
the multi-qudit case, albeit expressed here explicitly in 
terms of symplectic operators C (rather than invertible 
operators, which is only equivalent in the case n = 1). 

Corollary 5a: The action of U € C n (d) by conjugation 
on Pauli operators may be characterized by Eq. (|13l) . 

Proof: This is easy to show for Pauli operators, and 
follows for special Clifford operators from the character- 
ization in terms of Sp 2 „(Z/j). Then T;y is a product of 
linear operators arising from these special cases. □ 

While Refs. also relate the Clifford group to sym- 

plectic operators, we have shown more specifically that 
aC n (d) C C n (d) is a group representation of Sp 2n (Z£>). 
As well, because stabilizer tableaus transform linearly un- 
der Clifford operators (with the action on the Weyl block 
of the tableau being symplectic), we obtain: 

Corollary 5b: Acting on a proper tableau Tg by any 
conjugation tableau 7u yields another proper tableau. 



2. Generators for the symplectic Clifford group 

The symplectic Clifford group not only has convenient 
algebraic properties, but includes the best-known Clifford 
operators on qudits from the folklore, and analogues of 
operators from gen eralizations of the Pauli group to finite 
fields {e.g. Ref. 23]), including the single-qudit operators 



d-l 



d-l 



S = J2r q2 \q)(q\, F = ±= £ \p)(q\ , 



d-l 



(15a) 



M a = E NX<2l . for a e z *d > 

q=0 



where 7L* d stands for the multiplicative group of units 
modulo d and where multiplication is performed mod d\ 
as well as both of the two-qudit operators 



d-l d-l 



AZ 



E E r29192 i9i>M®i»>tei 



qi— qi — 



d-l d-l 



AX = E E ® \Q2 + 9i)(<Z2| 

Qi— qi — 



(15b) 



(15c) 



In the binary case, S is the 7r/4-phase gate and F is the 
Hadamard gate while AZ and AX are the controlled- Z 
and CNOT gates. For d ^ 2 in general, F is the quantum 
Fourier transform; and the operators M a , A, and AX can 
be used to perform any invertible affine transformation 
modulo d on the computational basis. 

For U unitary, we write A i > B to denote the rela- 
tion UA U'=B for the sake of brevity. One may easily 
verify that the above operators U are symplectic Clifford 
operators by computing the effect of UW&XJ^ on stan- 
dard basis states, for 1 < j < 2n, to verify the following 
equations: 





Z 








A 




l • 




Z 




1 




A 








z 




_1 (mod D) 

i 




A 


M^ x a 




Zd 


5 1 




5 1, 


Xi 


g) 1 


^A€ 


)Z, 


1 5 


■)Z 






1 5 


IX 




A, 


Zi 


5 1 


^Z£ 


5 1, 


A< 


E) 1 


^>A€ 


5 A, 


1 5 


5 Z 


AX^t 


®z, 


1 5 


DA 


^>1® 


A, 



(16a) 



(16b) 



(16c) 



all of which are mappings from Weyl operators to Weyl 
operators [23 |. Thus, each such U may be represented by 
conjugation tableaus of the form 1©Cj/, for the following 
symplectic operators Cjj € Sp 2n (Z/j) : 



L 1 

] 1 

1 

-1 



Caz = 



Cm — 



AX 



l l 

110 

10 

1. 

1-1 

10 

10 

11 



(17) 



Not only do these operators represent a collection of well- 
known operators which happen to belong to the symplec- 
tic Clifford group, they also characterize it: 

Lemma 6: The operators S, F, AZ (or AX), and M a for 
a ranging over the multiplicative units of Zc, generate 
aC n (d) up to global phase factors. 
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Proof: We show this result by showing that the operators 
CaXi C-s, Cp, and Cm o) together with direct sums with 
I2 correspondin g to tensor products with the identity op- 
erator on qudits [251, generate the group Sp 2n (Z£>). This 
result is technical, and we prove it by reduction to a simi- 
lar result shown in Ref. [H, Sec. IV] (see note [Hj])- Using 
Eq. (I16|) . one may show that AX = (1 ® Ft)AZ(l ® F), 
so that selecting either AX or AZ as a generator is equiv- 
alent. We may also describe AX2.1, a reversed ver- 
sion of AX 12 ■= AX (i.e. where AX 2j i has the second 
qudit the control, and the first qudit the target), by 
AX 2 ,i = (F f ® t)AZ(F (g) 1). It is easy to show that 
for g G Z, the operator AX^f has a conjugation tableau 
of the form 1 © C~ 9 , where 

"l " 



1 



measurement, obtaining a record h G Id of the eigenvalue 
r 2h associated to that eigenspace. 

For reasons that will become apparent in Section [TVl 
describing the evolution of a state under a non- 
destructive measurement is more complicated for d ^ 2 
arbitrary than for d prime. This is an obstacle to simulat- 
ing circuits which involve both multiple measurements, 
and unitary operations which are conditioned on mea- 
surement outcomes. We may overcome this by apply- 
ing the principle of deferred measurement, postponing 
all measurements to the end without affecting the mea- 
surement statistics (or final residual quantum states). 
This will allow us to efficiently simulate stabilizer cir- 
cuits with classically-controlled Pauli operations (but not 
classically-controlled Clifford operations in general) us- 
ing the techniques described thus far, so that we may 
consider the statistics of any single measurement uncon- 
ditionally from any measurement which precedes it. We 
may perform such simulations as follows. 



We can also generate SWAPi j2 gates on pairs of qudits, 
which interchange Z\ and Z 2 , and similarly X\ and X 2 : 
routine calculation will show that we may decompose 
SWAPi, 2 = (F 2 (g) l)AXi, 2 AX 21 AXi j2 • By construction, 
it has a conjugation map of the form 1 © for 



Cc 



10 
10 
1 
10 



(19) 



1. Reduction to Z r measurements 

As in the binary case [ll[ , there is no loss of generality 
in restricting from general Pauli measurements to "non- 
destructive" Z measurements (i.e. in which the measured 
qudit has a defined post-measurement state). For a W p G 
"P® n measurement, this may be done in a straightforward 
way by constructing a AW P ("controlled W p ") operator 
which on standard basis states \c) G Wd and |t) G 
would perform the operation 



Finally, the conjugation map of the operator FS~ 9 F> for 
g G Z is the operator 1 © CpCg 9 C~^ x , where one may 
compute 



1 
9 1 



(20) 



AW P \c) |t) 



w; |t> 



(21) 



We may construct W p G Vd on a single qudit straightfor- 
wardly out of AZ and AX gates. To illustrate this, con- 
sider the decomposition of W p into Z, X, and a "global 
phase" gate, as illustrated in Fig. [TJ 



The three operators of Eqs. (fTg)) - (j2"0|) . together with 
the operators Cp and Cm t , are precisely those used in 
Ref. @, Sec. IV] to decompose arbitrary symplectic op- 
erators modulo D. 

Thus an arbitrary conjugation tableau for an operator 
U G <sC n (d) may be generated by by those of the opera- 
tors S, F, AZ, and M a ; as U is characterized up to scalar 
factors by its conjugation tableau, the Lemma follows. □ 



C. Deferring measurements in stabilizer circuits 

For qudits of dimension d ^ 2, most Pauli operators 
are non-Hermitian. For a Pauli operator P of order at 
most d, we use the phrase "measurement of P " as short- 
hand, to refer to measurement of any Hermitian operator 
H with a spectral diameter less than 2ir, such that P = 
exp(iH). This amounts to collapsing the state of the 
system onto one of the eigenspaces of P via projective 



- W P1,P2 — = - ZP1 - T--P1P21 



FIG. 1. A notional circuit for a W p gate on one qudit. 



To obtain a coherently controlled version, given that 
W° = W cp = t-c'p^z^X ^, we replace the Z and 
X gates with AZ and AX gates on a common con- 
trol qudit, and act on the control with S~ PlP2 satisfy- 
ing S~ P1P2 |c) = t~ c P1P2 |c). Then we may decompose 
AW P = (S- piP2 © 1)AZP! AX^ 2 for single-qudit Weyl op- 
erators illustrated in Fig. [2] For a multi-qudit Weyl op- 



erator W n 



l,Pn+l 



Wp 2 ,p„ +2 



W Prl ^ P2rl where 



Pi ,p™+j 



act- 



n > 1 , we decompose AW P into operators AW 
ing with a common control and distinct targets. Using 
the Fourier transform F on the control qudit, we may 
apply standard techniques for eigenvalue estimation [271 ] 
to perform a W p measurement using a Z observable mea- 
surement on the control, as illustrated in Figure [3J 
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— X V2 - ZP* 



5-P1P2 — 



w u 



Wp 2 ,P„+2 



FIG. 2. Circuits for AW P gates, for W p £ V d and W p £ Vf 



|0) 



£iHD- w 



FIG. 3. Circuit to perform a measurement of an n-qudit Weyl 
operator W p , using a Z measurement on the top qudit to 
obtain the result. (This is illustrated here for measurement 
of a r 2h -eigenstate of W p , where the outcome will be h.) 



Using these restrictions, we may commute Z measure- 
ments on measurement registers past any Pauli operators 
which are conditioned on the outcome of that measure- 
ment. If a measurement record acts only as a classical 
control for Pauli operations, both the distribution of out- 
comes for r and the final state of any system depending 
on r is unaffected by commuting the Z r measurement 
past the controlled operation. More generally, while a 
controlled transformation depending on a measurement 
record s will transform Z r non-trivially, the group of ob- 
servables (Z m , . . . , Z r , Z s ) on the measurement records 
is preserved; so measuring Z m , . . . , Z s after the trans- 
formations will yield the same distribution of outcomes, 
as measuring Z m , . . . , Z s first and then transforming the 
outcomes. With either approach, one may simulate all 
classically controlled operations P € p® n which depend 
on measurement outcomes using the corresponding co- 
herently controlled operators AP. The distribution of 
any single measurement outcome may then be simulated 
using the techniques described thus far. 



Remarks on classically controlled Clifford operations 



Measurement of a Pauli operator P — t W p may be 
performed similarly by replacing AW P with a controlled- 
P gate, AP = [Z 5 ® 1®") AW P . These substitutions may 
be easily performed at the time of simulation. 

2. Qudit states as measurement records 

While it is common to conceive of single-qudit mea- 
surements as being destructive, replacing the measured 
qudit with a classical record, we consider the qudit r on 
which we perform the Z observable measurement to be a 
persistent quantum system whose state is the measure- 
ment record. To retain consistency with the notion of a 
classical measurement record, we may impose constraints 
of the following sort: 

• We may require that no operators act on the qudit 
r (after the measurement) aside from those which 
commute with Z r . For instance, it may act as the 
control of a product of AX or AZ operators, instead 
of acting as a classical control for Pauli operators 
on the same target qudits. In this case, the state 
of the qudit r remains an unchanging record of the 
measurement outcome. 

• We may instead allow measured qudits r to be 
subject to Pauli operators themselves, with clas- 
sical control limited only to qudits s which are also 
measurement records or input parameters. Such 
transformations of the measured qubits are limited 
to invertible transformations of the standard ba- 
sis over Zrf, and may be simulated in the stabilizer 
formalism by representing them as Clifford group 
operators. 



Deferring measurements cannot be done for measure- 
ment records which control Clifford operations while re- 
maining in the stabilizer formalism, using the techniques 
described above, except for controlled-Pauli operations. 
Even in the simplest case, acting on states in which the 
control is in a standard basis state, such a simulation 
would involve transforming some columns of the tableau 
(i.e. representing Pauli operators supported on the target 
qudits of the operation) with different linear transforma- 
tions depending on the values of other columns (which 
represent Pauli operators supported on the control qudit 
of the operation). Such a transformation is at best multi- 
linear, but not linear. Furthermore, this transformation 
would be ill-defined for any column in which the control 
qudit is acted on by an operator with a non-zero X com- 
ponent. For instance, the coherently controlled operator 
AS obtained by deferring the measurement on the control 
acts on the operator Iglas follows: 

AS(X <g> l)AS f |c) \t) 

\ (X <8> S* 1 ^) |c) \t) , ifc = <f-l. 

Thus, AS fails to map X ® 1 to another Pauli operator 
via conjugation, and isn't represented by a transforma- 
tion of stabilizer tableaus. (Indeed, as it fails to preserve 
the Clifford group, such an operator AS is not itself a 
Clifford operator; supplementing the Clifford group with 
any single non-Clifford unitary generates a set of oper- 
ators approximately universal for quantum computation 
for d — 2 [28| as well as any other prime dimension [29|, 
Appendix D], so that classical algorithms to simulate 
toffoli gates together with Clifford operations are im- 
possible unless BPP = BQP.) To simulate a Clifford 



8 



operation controlled by measurement outcomes, we must 
therefore explicitly compute (the probability distribution 
of) the measurement outcome, and use this outcome to 
determine mid-simulation which linear transformation to 
apply to the target qudits. 

For any fixed measurement outcome, or for a classi- 
cal control value provided as input (or generated accord- 
ing to a fixed probability distribution) , simulating a clas- 
sically controlled Clifford operation is no more difficult 
than a fixed Clifford operation; the only obstacle in this 
case is simulating the evolution of the state caused by 
the measurement itself, which we treat in Section IfVI 



D. Computing the distribution of a single terminal 
measurement 

Describing the evolution of a state under measurement 
proves to be more technically complicated in the case of 
composite dimensions d ^ 2, as we shall see in Section llVl 
However, determining the distribution of outcomes for a 
single measurement, on a state which can be represented 
by a proper stabilizer tableau, can be treated using the 
techniques described thus far. We now summarize how 
to determine the distribution over possible outcomes, de- 
ferring the account of how states transform under mea- 
surement for Section HVl and proofs of certain details to 
Appendix [B] 

Consider a proper stabilizer tableau Tg which repre- 
sents a group Gg stabilizing a unique state. To simulate 
measurement with respect to a Pauli operator P, con- 
sider some Weyl operator W p cx P. While W p may not 
commute with all elements of Gg, there may be a non- 
trivial power of W p which does. The smallest positive 
power s > f such that W p does commute with all of 
Gg constrains the distribution of possible measurement 
outcomes for P. By Lemma [TJ s is the smallest positive 
integer such that [sp, v&] = (mod d), where ipk © vj, is 
the k th column of T§. That is, s is the smallest positive 
integer such that 

s(0®p) T (0® ct 2 „)T s = T (modrf). (23) 

From the minimality of s, it follows that s will be a factor 
of d, with s = 1 in the case that W p commutes with Gg. 
Specifically, let (f> = (0 p) T (0 © <J2n) Tg, and consider 
the greatest common factor rj of the coefficients of <p to- 
gether with d: then s — d/rj. As commutes with all 
of the generators of Gg, and as Gg stabilizes a unique 
state, we may show (Lemma H2]) that W p is proportional 
to an clement of Gg. Then the pre- measurement state 
is an eigenstate of W p , whose eigenvalue constrains the 
possible measurement outcomes. If it is a r 2 *-eigenstate, 
then r~ 2t W sp stabilizes the state, in which case there 
is a vector t' sp' in the column-span of Tg such that 
T~ 2t W sp — r~ 2t W sp >. In order for this equality to 
hold, we require that sp' — sp = (mod d), and that 
2t' -2t = [sp', sp] (mod D) as well. 



• For d odd, these two conditions imply that t' — t = 
(mod d), so that we simply need to determine for 
which value of t £ that t © sp is in the column 
span of Tg . 

• For d even, the vector i'©sp' may differ from t©sp. 
In particular, we may have t' and t may differ by | 
modulo d; and the vector dx = sp' — sp corresponds 
to some representation of the identity operator 1 = 
Wdx which accounts for the difference in the phase 
coefficients, by the formula 

T- 2t 'W sp , = T- 2t +[ dX < S P] W sp+dx 

= r- 2t W sp W dx , (24) 

by Lemma [2j We may attempt to account for the 
contribution of Wd-x. to the phases by introducing 
auxiliary columns of the form Uj := 5 [ej, sp] © dej 
to the tableau Tg : adding Uj to (t © p) yields an- 
other Pauli vector representing the same operator. 
We may also introduce a column uo = (d © 0), 
also representing the identity operator, to account 
for differences of ±d in the phase coefficient which 
are insignificant. Thus, it suffices to determine for 
which value of t £ 7Ld that (t © sp) is in the column 
span of [u ui • • • u 2 „ | Tg ] . 

In either case, it suffices to adjoin the vector — (Offisp) to 
the tableau (together with the vectors Uj in the case of d 
even), and perform column-reductions to clear the Weyl 
tableau. The resulting matrix will contain a vector (t®0), 
which implies that (t sp) is the the column span of the 
original matrix and that the pre-measurement state is a 
r 2t -eigenstate of W p . 

The constraint that this imposes on the outcomes of a 
W p measurement are as follows. As (r _2t I / Fp) r) is pro- 
portional to 1 and has +l-eigenstates, we have 2t-q = 
(mod D), from which it follows that t is a multiple of 
s. The post-measurement state is a r 2tl -eigenstate of W p 
by definition for some u £ 2j, and in particular must 
be a r 2 *-eigenstate of W p as noted above; thus 2t = 2su 
(mod -D), or equivalently u = t/s (mod rj). We show in 
Appendix lB 3l that the outcomes are in fact uniformly dis- 
tributed over all u £ Id satisfying this constraint. This 
characterizes the probability distribution over the out- 
come r of the measurement of W p . For P = t~ 2S W p , 
the distribution of outcomes for a P measurement are 
uniform over the solutions h £ 7*d to the congruence 

h = u + S = trj/d + S (mod 77), (25) 

where again 77 = gcd(d, <f>i, (f>%, . . .). In particular, the 
results are uniformly distributed over some coset K + rjZd, 
where we may take k = tr]/d+ 5. 

Note that the above analysis accommodates the pos- 
sibility of a deterministic outcome. If P commutes with 
Gg, then r\ — d, in which case Eq. (|25|) has a unique so- 
lution mod d: the outcomes are "uniformly" distributed 
over a singleton set (i.e. follow a delta-distribution). This 
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analysis also generalizes the familiar scenario of a random 
measurement outcome for d prime: in that special case, 
we have r\ = 1 if the measurement observable P does 
not commute with all elements of S. In the case r\ = 1, 
Eq. (|25p trivializes and imposes no constraints on h G Z^, 
and so is uniformly distributed over Z^. 



E. Summary for simulating unitary stabilizer 
circuits 



We have shown that using Weyl operators to describe 
stabilizer tableaus allows us to describe a simple, linear 
formalism for representing the action of Clifford circuits 
on standard basis states, and determining the measure- 
ment outcomes for single measurements. These results 
follow from the special algebraic properties of Weyl op- 
erators as described in Lemma [TJ and the fact that stan- 
dard basis states may be represented by proper tableaus 
(as defined in Section fill Al) . Demonstrating a subgroup 
aC n (d) C C n (d) which is a representation of Sp 2n (Z£>) is 
an additional dividend of this formalism. 

The remainder of the article is devoted to extending 
this formalism to explicitly describe the evolution of sta- 
bilizer states under measurements, which involves provid- 
ing a linear analysis of those states which lack "proper" 
stabilizer tableaus as described in Section IIII Al 



and inconsequential changes in the phase block): 



2 
2 



2 
6 



6 
2 



6 
6 



(26) 



where each * may independently be either or 4 (rep- 
resenting a global phase of +1). The reader may verify 
that none of these are proper, as the columns of the Weyl 
blocks are not orthogonal under the symplectic product 
modulo D — 8. Thus, there do not exist any proper sta- 
bilizer tableaus for §; by linearity, there are no proper 
tableaus for any other set of operators generating the 
same group either. 

As the Pauli vectors of commuting operators are al- 
ways orthogonal mod d under the symplectic inner prod- 
uct, such obstacles do not arise for d odd. In the case of 
d = 2, while "improper" tableaus do exist, every stabi- 
lizer state has a proper tableau, which may be obtained 
from an improper one by the techniques of Appendix IA II 
However, to describe a formalism which functions for ar- 
bitrary d, it is necessary to incorporate corrections to the 
phases as in the standard "binary" formalism, in order to 
describe different generating sets of the stabilizer group. 
We show that these corrections may be subsumed in a 
more general linear formalism, as follows. 



Extended stabilizer tableaus 



IV. IMPROPER TABLEAUS AND EVOLUTION 
UNDER MEASUREMENTS 

As we note on page 2] (and show below), there do not 
always exist proper stabilizer tableaus for a stabilizer 
state in the case of even d > 2. Furthermore, measure- 
ments may transform a state stabilized by a group with 
a proper tableau, to one stabilized by a group without 
any proper tableau. This presents an obstacle for a uni- 
form treatment of stabilizer circuits in the manner we 
have described so far, and are unavoidable when simu- 
lating the distributions of outcomes of circuits involving 
multiple measurements for even d > 2. We must there- 
fore extend the formalism described thus far to describe a 
linear formalism to simulate arbitrary stabilizer circuits, 
and simulate outcomes of multiple measurements. 



A. Extension to accommodate states having no 
proper tableaus 

For even dimensions d > 2, there are stabilizer states 
whose stabilizer groups cannot be represented by a 
proper stabilizer tableau. The simplest example is the 
single-qudit state ^j(|0) + |2)) for d = 4, which is sta- 
bilized by the set § = { Z 2 . X 2 } . This set of generators 
has four possible stabilizer tableaus (up to column swaps 



Consider a set of commuting Pauli vectors p 1; . . . , p^, 



where each vector p 



represents a Pauli op- 



erator Pj = t 'W Vj . By Lemma [TJ we have W U W V = 
7-[ u ' v lW u + v , where [u,v] = [v, u] as Ph and Pj com- 
mute. It follows that [u,v] € {0, d} modulo D. These 
represent corrections to the phase coefficients of Pauli 
vectors by multiples of ^ (which would always be in the 
case of d odd). When computing products of the Pauli 
operators Pj, we may avoid computing the quadratic 
dependency arising from the symplectic inner products 
[v/i, Vj] by storing an array Sg {0, |, d, ^ } ixi such that 
Ehj = \ [v/,,Vj] (mod D). We may do this using an 
expanded table of coefficients, of the form 



T s = 



vi 



Vf 



(27) 



We call such an array an extended stabilizer tableau, and 
refer to 3 as the phase correction block of the tableau; 
the tableau T§ is proper if 3 = (mod D). In the case 
of d odd, this in fact holds by necessity. 

Products and recombinations of generators may be 
represented by linear super operators on such extended 
tableaus, differing only slightly from left- and right- 
matrix- multiplications on T§. For each 1 ^ j ^ £, let 
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Pj = Pj © Sej be the column of T§ extending . Con- 
sider a product of Pauli operators: 



P 



PhP, 



-2(4> h + 4> j +B. Kj ) j^r 



(28) 



Let tp = cj>h + 4>j+ ^h.j for the sake of brevity. When tak- 
ing products of P with other operators , the quadratic 
dependency involved in the phases depends only on the 
Weyl operator W Vh + Vj , to which P is proportional. 
Therefore we have 



pp. — T--2(v + <Ak + [v h +v J ,v i .]/2) tit 



2(if + 4> k + [E h}k + 5j,k]) 



(V h +Vj)+V fc ) 



(29) 



where the phase correction [S^fe + z,j,k] is a sum of the 
phase corrections arising from Ph and Pj . Then, we may 
represent operator P by an extended column vector of 
the form 



-h,j 



Ph + Pj 



-h,j 



(30) 



Thus to obtain the extended tableau Tg in which the 
generator Pj is replaced with P = PhPj, it suffices to 
add the h th column into the j th column, subtract the h th 
row (of the phase-correction block) from the j th row to 
maintain the antisymmetry of that submatrix, and finally 
add the corrective term Shj into the phase coefficient. 
That is, we compute 



T 5 = * 



h-J 



E. 



2n+l+h,2n+l+j Eh,j 



(31a) 



where E a ,b = 1 + e/,e a ^ s the usual elementary matrix 
for adding row a into row b, and where 



*hj(M) = M 



eie 



2n+l+h 



M 



(31b) 



is a superoperator which adds the appropriate phase cor- 
rection (stored in the j th column) for the operators rep- 
resented by columns h and j, to the phase coefficient of 
column j. 



2. Reducing the extended tableau modulo d 

In the case that d is odd, we have D — d, and the entire 
tableau remains reduced modulo d at all times. However, 
if d is even, coefficients c € {d, d+1, . . . , 2c?— 1} may arise. 
Though these coefficients never present any actual diffi- 
culties in representing stabilizer states, it is still possible 
to reduce such coefficients modulo d by performing phase 
corrections if this is desired. As an extended tableau is 
under no restrictions to remain proper, we may more 
freely choose which Pauli vectors we use to represent a 
given Pauli operator generating a stabilizer group, so long 
as we maintain the correct phase correction information. 



For a single Pauli vector, we may reduce the Weyl block 
modulo d by considering different representations of the 
identity operator 1 by Pauli vectors, such as 



W d&1 
W d * 



xf. 



= x° 



(32) 



This corresponds to Pauli vector representations of the 
form © d&h for 1 $5 h ^ 2n. The phase correction 
required to combine such a vector with an arbitrary Pauli 
vector tp © v is given by 



2 V h+n , 



-,Vh-r 



if 1 ^ h ^ n; 
if n < h < In, 



(33) 



modulo D. Because — | = d+ 1 (mod D) for d even, and 
differences of d have no effect in the phase columns, the 
phase corrections of — s may be replaced by phase cor- 
rections of +5, or vice versa. Reduction modulo d of the 
Weyl coefficients of Pauli vectors can then be performed 
with the transformations 



ip®v i — ► ((p+^v n+h )®(v-de h ), 
ip®v i — > {ip + iv h ) © (v - de n+h ), 



(34) 



as the vectors on the left and right all represent the same 
Pauli operator t~ 2v W v . (This is used implicitly in the 
algorithm for determining distributions of measurement 
outcomes in Section fill Dl ) By subsequently reducing the 
phase coefficient modulo d, the entire vector may be re- 
duced modulo d. 

While the above transformations do preserve which 
Pauli operator is represented by a given Pauli vector, 
it does not preserve the symplectic inner product of the 
Weyl block with those of other Pauli vectors. To main- 
tain an extended tableau while reducing Weyl coefficients 
modulo d, one must also update the phase correction 
block S. Suppose that ipj © Vj © Se,- is the j th column 
of an extended tableau Tg which represents a set of op- 
erators {Si, . . . , Sj, . . . , Si}. As with the operation to 
represent multiplication of generators described in Sec- 
tion IIV A 11 we may reduce the h th coefficient of v by 
d by combining the phase correction information of the 
vector © det for the other generators Si, ... ,St, to the 
vector £ describing the phase correction information of 
(p ffi v. That is, we compute the row vector 



$h = 



[de h ,vi] ■■■ [de h ,v e ] 



\eJ l [o\a 2n \0]T s (mod d), 



(35) 



and then add £ h to the j row of the phase correction 
block S (updating the phase correction information for 
all other columns with respect to the j column), and 
subtract ^ from £ (updating the phase correction in- 
formation for the j th column itself). Note that by sub- 
tracting £^ from £, the phase correction necessary for 
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v is added into the j th coefficient of £; before updat- 
ing the j th row of the entire phase correction block, we 
may apply the superoperator iffjj to perform the neces- 
sary phase correction in the j th column. As the coeffi- 
cients of the phase correction block are only ever added 
to the phase coefficients, they may themselves always be 
reduced modulo d as in Eq. (|35[) ; thus whether the vector 
£ h and its transpose are added or subtracted is immate- 
rial, except to maintain the antisymmetric property of 
the phase correction block. Thus, to reduce the (h,j) 
coefficient of the Weyl block of an extended stabilizer 



tableau for 1 ^ h ^ 
compute 



and obtain a new tableau X§, we 



-3 



(0 © © H)eJ - de 1+h e] 



— ( '' ; .. (^s + ^e 2n +i+jeJ l+1+h T s ^ 

- |(0 © © [T s T e„ +1+ ,,.])eJ - de 1+h ej ; (36a) 

to reduce the (h+n,j) coefficients of the Weyl block for 
1 ^ h ^ n, we instead compute 



&2n+l+j$,h 

(0 © © £[)e] - de n+1 
|(0ffi 0© [Tje 1+h ])ej -de. 



-l+h&j . (36b) 



In both cases, the mapping is an affine superoperator, 
involving both the superoperator * and transposition of 
the matrix T§, and also of course the translation opera- 
tion by — dei+h&J or — de n +i+h&j to reduce the actual 
coefficient of the tableau. 



3. Fixing the number of rows of extended stabilizer tableaus 

When simulating unitary transformations on states 
represented by extended stabilizer tableaus, it may be 
preferable to use linear operators of a fixed size, inde- 
pendently of the state being transformed. While the 
phase correction block is a square symmetric matrix, and 
may have a variable number of rows, we may bound its 
size if we can establish a maximum number of generators 
needed to describe an arbitrary stabilizer group. We can 
then pad the number of rows or columns in any tableau 
on n qudits to obtain tableaus with a fixed size. 

As noted in Refs. [HEB], a stabilizer group on n qudits 
may have a minimal generating set strictly larger than n 
when d is composite. The simplest example is the single- 
qudit state — ^ ( 1 0) + |2)) described above for d = 4, for 

which {Z 2 ,X 2 } is a minimal set of generators. (Similar 
examples exist for any composite d > 2; one may consider 
a state stabilized by {Z dl , X d2 } for arbitrary di,d 2 > 1 



on qudits of dimension d — d\d 2 .) However, we may show 
that a minimal generating set for any n-qudit stabilizer 
group has at most 2n generators, as follows. For any 
extended stabilizer tableau T§ (proper or otherwise) with 
I > 2n columns, the Smith normal form Sg = L W§i? of 
its Weyl block W§ has the form 



h 
t 2 











tan-i 
t 2n 







(37) 



where each coefficient tj is either zero or a divisor of 
D, and where tj + \ is a multiple of tj for each j < 2n. 
The right- most (£ — 2n) columns of W§i? are then also 
zero. Interpreting R as a recombination of the columns 
(and applying the appropriate phase corrections via 
Eq. pial) - ). we obtain an extended tableau Ts which has 
W§i? as a Weyl block. The right-most (£ — 2n) columns 
then each represent the identity operator, and may be 
omitted to leave a tableau with at most 2n columns. As 
the phase correction block requires only as many rows 
as columns, we may specify an arbitrary stabilizer group 
with an extended stabilizer tableau having at most 4n+l 
rows: one for the phase block, and 2n each for the Weyl 
and phase-correction blocks. Thus, by padding, we may 
fix extended tableaus for n qudit states to have An + 1 
rows and 2n columns. 



Simulating unitary stabilizer circuits on 
extended tableaus 



We may extend conjugation tableaus, representing the 
effect of conjugation of Pauli operators by Clifford op- 
erators, very simply. As symplectic transformations do 
not affect the pair-wise symplectic inner products of the 
vectors Vj, we may leave the phase correction blocks of 
the vectors untouched when we describe a transfor- 
mation by U € C n (d). An extended conjugation tableau 
corresponding to a conjugation tableau Tjj can thus be 
given by a square matrix 



'7u 








l2n 



' 1 


hi ■ ■ ■ h2„ 


••• " 











Cu 










l2n 


. 







(38) 



where again Cjj is a symplectic transformation of Zj^ . 

Simulating terminal measurements as in Section [ill Dl 
requires only the ability to compute alternative generat- 
ing sets for the stabilizer group as described in Eq. (|31aD , 
and the ability to compute symplectic inner products as 
in Eq. (j2U)) . For the latter, symplectic inner products of 
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pairs of operators v/,, Vj are not pertinent for computing 
the inner products [p, v h ] and [p, v 3 -] which are necessary 
to determine the statistics for measuring an observable 
P = T~ 2v W p . Thus it suffices to compute the smallest 
integer s > 1 such that 



sp 











T s = T (mod d) 



(39) 



as a simple modification of Eq. (|23[) . We may then carry 
out a simple variation of analysis of Section !!!! Dl defining 
£h as in Eq. ([33]) , we extend the tableau Tg with vectors 
iih = © deh © £j for 1 j ; $C 2n in place of the vectors 
u/j and u' h described in Section IIII Dl and attempt to 
clear the Weyl block of the vector — (0 © sp © 0), using 
Eq. (|31a[) to perform column combinations. 

Adopting these extensions to stabilizer tableaus and 
conjugation tableaus allow us to simulate arbitrary sta- 
bilizer states, whose stabilizer groups may or may not be 
representable by a proper tableau, using linear operations 
over Zjj as in Section iHTl 



Simulating stabilizer state evolution under 
Pauli measurements 



In many cases, including the case d = 2 and d odd, 
every stabilizer state has a proper tableau, and simulat- 
ing Pauli measurements in general may be achieved by 
transformations of proper tableaus. However, for d > 2 
even, such measurements may transform a state from one 
which can be represented by a proper stabilizer tableau, 
to one which cannot. For instance, in the case d — 4, we 
may prepare the state 

NA) = AX a 2 b F a |0)jO) b 

= U\°) + \ 2 ))a®\°) b + MW+IS))^)*,, (40) 

for which we may easily construct a proper tableau us- 
ing the techniques of Section 11111 However, if we measure 
Zf, on \ip), the post- measurement state on a would differ 
from (|0) + 1 2)) by at most an X operation. As we 
saw in Section IIV Al such states have no proper stabi- 
lizer tableaus. We must therefore extend the formalism 
of Section Hill to deal with evolution of states under mea- 
surement in the general case. 

We now describe a procedure to describe this evolution 
for all d 2, using the extended tableaus introduced in 
Section TlV Al deferring the proofs of certain statements 
to Appendix [B] (In Section IIV C 31 we describe how this 
analysis may be simplified in the cases d = 2 and d odd, 
for which extended tableaus are not in fact needed.) Us- 
ing the constructions of Sections IIII CI and IIVB1 we may 
reduce the problem to that of describing the effect of a 
single-qudit Z measurement. (Again, this has the side- 
effect of allowing us to simulate classically controlled 
Pauli operations depending on the measurement result 
using a coherent controlled-Pauli operator. We sketch a 
more direct approach to measurement in Section [IV C 31 ) 



Let S = {Si, . . . , Si} be a minimal generating set 
for a stabilizer group Gg which has a unique joint +1- 
eigenstate. A Z r measurement has a non-trivial effect 
on a stabilizer state, and a non-deterministic outcome, if 
and only if the state is stabilized by a Pauli operator that 
does not commute with Z r . Given an extended stabilizer 
tableau T§ , we may compute whether Z r commutes with 
each operator Sk £ S by computing 



= K 











T§ = (0 © e n+r © 0) T T S (41) 



following the description in Eq. (|39[) . The result is a row- 
vector 4> = [0i ... <f> £ ]e Z£ such that Z r S k = T 2<t, *S k Z r . 
As the row- vector <j> is the (n + r+ l) st row of the tableau 
T§, the tableau has the following block structure: 



row 1 
row 2 





Ti,i ■ ■ ■ Tij 




T 2 ,i ■ ■ ■ T 2l £ 


Ts = 


01 • ■ ■ 4>f 











(42) 



■ row n + r + 1 
row In + 1 

phase correction block 



We may recombine the generators described by the 
columns using Eq. (j31a[) . to obtain a tableau Tg which 
has at most one column — without loss of generality, the 
£ th column — containing an entry r] ^ (mod d) in the 
(n + r + l) st row: 





T M ■ 




T14 


4— row 1 




T 2 ,i ■ 


T24-1 


T 2 ,e 


<— row 2 


Ts = 








V 


4— row n + r + 1 




72n+l,l ■ 


■ ^Wf-i,/— 1 




A— row 2n + 1 









(43) 



We may restrict rj to be a divisor of d if it is non-zero, by 
multiplying the £ th column by a suitable scalar a £ Z* D . 
(To maintain the antisymmetry of the phase correction 
block, the I th row of S must also be multiplied by a 
in this case.) Note that rj £ {0,1,2} in the case that 
d — 2, T) £ {0, 1} in the case that d is an odd prime, and 
rj £ {0, d} for any d ^ 2 in the case that Z r commutes 
with every generator of S. 

The tableau Tg describes generators S = {Si, . . . , Si} 
for Gg in which only Si may fail to commute with the 
observable Z r , depending on the precise value of 77. Let 
Vj represent the Weyl block of the j th column of Tg , and 
let s = d/rj (if rj > 0), or s = 1 (if rj = 0). Then s > 1 
is the smallest power of Si such that Sf commutes with 
Z r , as [sv£,e r ] = (mod d) by construction. Similarly, 
Z^ commutes with Se, and thus with every element of S. 
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1. Computing the distribution of measurement outcomes 

As G§ stabilizes a unique state, it follows (by 
Lemma IT2"j) that Z^ is proportional to some element of 
G§ . We introduce the vector 

i ,t [ — , if d is even and ri > 0: ,,,, 
10, if a is odd or 77 = 0, 

for the sake of brevity. This vector then contains the 
phase correction information needed in order to represent 
products of Z^, which we represent by the Pauli vector 
(0 © se r ), with the generators Sj in an extended tableau. 
We find the element of Gs to which Z^ is proportional 
using a similar technique to Section IIVB1 temporarily 
adjoining columns = © de.h ® $J f° r £/i as given in 
Eq. (|35p to aid in performing column reduction on the 
tableau, using Eq. (|31a[) to perform column recombina- 
tions, to clear the Weyl block of the vector — (0ffise r ©£). 
The result will be some vector t © © 0; by linearity, it 
follows that the operator 

R := r- 2t Z s r = T- 2t W S e r , (45) 

represented by the extended Pauli vector t © se r © £, is 
the operator which is proportional to Z^. and contained in 
Gs- As in Section llll D[ the coefficient t must be a multi- 
ple of s in order for R to have a non-trivial +l-eigenspace. 
Because R commutes with Gs and with Z r , it stabilizes 
both the pre-measurement and the post-measurement 
state. The operator R = T~ 2he+1 Z r which stabilizes r 
after measurement satisfies R s = T^ 2hs Z^ = R; that 
is, the measurement outcome h must satisfy he+i = t/s 
(mod 77), or he+i € t/s + r\Ls,. We may show (see Ap- 
pendix IB 3|) that the outcome is in fact uniformly dis- 
tributed among the residues modulo d satisfying this con- 
straint. 

In the case s = 1 (i.e. if 77 € {0,d}), the state is sta- 
bilized by some operator r~ 2h Z r prior to measurement; 
the outcome is then a delta-peaked distribution, or "uni- 
formly distributed" across the coset h + r/Zd = {/i} C Zj. 
Furthermore, the measurement of Z r does not affect 
the state of the system; as an alternative to the above, 
we may transform the tableau to explicitly represent 
the fact that some operator r~ 2h Z r stabilizes the pre- 
measurement state by column-reduction, using the op- 
erations for column combinations and reduction modulo 
d of Sections IIVA1I and IIV A 21 The resulting tableau 
will contain a column of the form h © e r © £ explicitly 
representing the operator r~ 2h Z r . 

2. Computing and transforming between post-measurement 

states 

Having fixed a given outcome he+i of the Z r measure- 
ment, we perform the following transformations on Tg 
to represent the post-measurement state. As the state is 



unchanged if Z r commutes with all stabilizer generators, 
we restrict ourselves to the case s > 1. 

As Se does not commute with Z r , it does not stabilize 
the post-measurement state. But by construction, Sf 
does. To represent this, we multiply the £ th column of 
T§ by s. To keep the phase correction block consistent, 
we multiply the £ th row of the phase correction block by 
s as well. We then augment the generating set to include 
the operator R stabilizing the measured qudit: we do so 
by replacing the adjoined (£ + l) st column of the tableau 
with /i£ + iffie r ffi£ . (Note that by construction, £ as given 
in Eq. (|4"4"| is also the appropriate phase correction vector 
in this case). We must also extend the phase-correction 
block by an additional row to represent the new generator 
of the stabilizer group; by symmetry, this is simply ^ej 
if d is even, or T for d odd. This yields a new extended 
tableau Ts< representing a set §' of stabilizer generators 
for the new state. 

It may occur that the operator Sf represented by the 
£ th column of f§> is in fact 1, having a Weyl block con- 
sisting only of multiples of d. This will occur for instance 
if 77 = 1 (which always occurs in the familiar case of d 
prime) . If this is the case, we may drop the £ th column 
of the new tableau entirely, and over- write the £ th row of 
the phase correction block with the (£ + l) st row rather 
than filling a new row. More generally, it may be the 
case that Sf can be generated by the other operators 
{S[, . . . , S , £_ l , S' i+1 } C §'; this necessarily occurs when 
£ = 2n, as this is the maximum number of operators 
needed to generate the group Gs' ■ To discover whether 
this is the case, we may attempt to clear the Weyl block 
of the £ th column of T§f modulo d (i.e. reduce it to a 
{0, <i}-vector), as always using column operations with 
phase corrections as in Eq. (|31ap . If successful, the re- 
sulting vector in the £ th column represents the identity, 
and may be discarded as described above. 

To simulate the transformation which occurs for a par- 
ticular measurement outcome h* € t/s + r/Ld-, we may 
simply set he+± := h* by fiat, and perform the opera- 
tions as above. In a physical setting with actual stabi- 
lizer circuits acting on qudits of dimension d, as in the 
familiar case d = 2, one may simulate fixing the outcome 
by transforming the post-measurement state unitarily be- 
tween the possible outcomes. For d 2 arbitrary, we may 
show that it suffices to act on the post-measurement state 
with the operator Si, which is represented by conjugat- 
ing each of the elements of §' by Se,. By construction, 
this operator commutes with the generators Sj = Sj for 

1 ^ j < £, as well as the generator S' e = Sf, which are 
represented by the first £ columns of the tableau T§/; it 
only fails to commute with the generator S[ +1 = R oc Z r , 
transforming it to R' = r~ 2, 'i? = r - 2 ('^+i+')) Z r instead. 
Repeated applications of powers of Se then suffice to pro- 
duce any desired post-measurement state by adding a 
suitable multiple of 77 to the exponent in the phase. 

The entirety of the above analysis, for non-commuting 
Z r measurements, generalizes the well-known stabilizer 



14 



formalism in prime dimension (as in Ref. [8|). If d is 
prime, we have 77 = 1 as we have noted above, so that the 
post-measurement state is no longer stabilized by St but 
instead by R = r~ 2h Z r , for some h uniformly distributed 
over Z,j. The novel features for d arbitrary are that the 
possible outcomes over which h varies can in general be 
any coset of the form k + r(Ld (which is equal to the 
whole group Id only when 77 = 1), and that the post- 
measurement state is still stabilized by (which is 
trivial only when rj = 1). 

3. Special cases permitting simplified evolution under 
measurements 

The extensions of the preceding sections to the sta- 
bilizer formalism allow us to simulate arbitrary stabi- 
lizer circuits for arbitrary dimensions d 2. We now 
briefly remark on simplifications which are possible in 
some special cases, which may lead to modest savings in 
the amount of effort and work-space required in practice. 

When simulating a measurement of a Pauli operator 
P = t~ 2S W p , it may be that no further operations de- 
pend on the measurement outcome (i.e. on the eigenvalue 
r 2h corresponding to the post-measurement state); or if 
the outcome is used only to control Clifford operations, 
that representing the outcome by the state of a qudit 
provides no advantage. In either case, the introduction 
of an explicit measurement register for controlled-Pauli 
operators to act upon is unnecessary. We may compute 
the distribution of outcomes, and the post-measurement 
state conditioned on any particular outcome, in a similar 
way as described above but without the introduction of a 
measurement register. Elaborating Eq. (|4ip . we compute 

= (0 8 p ® 0) T !T S , (46) 

and consider whether it is equivalent to zero modulo d; if 
not, we consider what reversible column-transformation 
operations would map it to a row- vector of the form rjej , 
and then apply those same transformations to the tableau 
Tg. Doing so yields a tableau representing a generating 
set S = {Si, . . . , Si} in which only Si fails to commute 
with P; one may easily generalize the analysis above from 
that point on, substituting the Pauli vector (0©e r ) repre- 
senting Z r with the vector (<5©p) representing P. Doing 
this is substantially similar to performing the same op- 
erations as in Section IIV C 21 as performed on a tableau 
where we have simulated the P measurement by a Z r 
measurement, with the primary difference being that we 
omit the additional column involved by explicitly intro- 
ducing the ancilla r. 

As noted in Section [HI Al extended stabilizer tableaus 
are unnecessary in the case of d either prime or odd, as 
all stabilizer groups may be represented in those cases 
by proper stabilizer tableaus. In the case of odd d, no 
special effort is necessary, as all stabilizer tableaus are 
proper in that case; we may omit phase correction blocks 
in that case. In the remaining case d = 2 (i.e., for qubits), 



further effort is required to ensure that the tableau of a 
post-measurement state is proper if we are to dispense 
with phase correction blocks. If at least one of the ele- 
ments of the stabilizer group Gg fails to commute with 
the Pauli operator P = t~ 2S W p being measured, we may 
as usual obtain a generating set {Si,..., St} in which 
exactly one generator St does not commute with P; it 
suffices to compute a vector p = p (mod 2) such that 
P = t~ 2S W p — T~ 2S Wp, which is orthogonal modulo D 
(where D = 4 in this case) to the columns of the Weyl 
block Wg representing Sj for 1 ^ j < i. We may com- 
pute p using the techniques presented in Appendix IA II 
in this case, which should have a single solution. If how- 
ever P commutes with the entire stabilizer group, there 
is no choice in how it may be represented in order to 
maintain a proper tableau: an operator proportional to 
P is already generated by the group, and the system of 
equations in Appendix IA II determining a suitable repre- 
sentative for that operator has a unique solution p, for 
which [p, p] = 2 (mod 4) may hold. In this case, how- 
ever, no transformation of the state occurs upon measure- 
ment, so the existing proper tableau suffices to describe 
the post-measurement state. We may then dispense with 
extended tableaus in these cases if desired (though ex- 
tended tableaus still provide the benefit of making pos- 
sible reduction of the coefficients modulo 2 for tableaus 
over qubits). 

Finally, for d prime, certain elements of the analysis 
in Section IIII Dl and Section IIV CI may be simplified to 
yield the known results for simulations of measurements 
in those cases @. For instance, in the case that not all 
generators Sj £ S commute with the measurement opera- 
tor P, it is not necessary to perform column transforma- 
tions to obtain a column which represents an operator 
Si such that [P, St] — t 2 '' = t 2 ■ We may instead find 
any single operator Sj which fails to commute with P, 
and compute some non-trivial power Sj of it such that 
[P, Sj] = t 2 , and use this to obtain a generating set in 
which only Sj fails to commute with P by combining 
other generators with appropriate powers of S"j. Also, 
no non-trivial power of Sj will be represented by the ta- 
ble T$i for the post-measurement state; in particular, the 
tableau will never have more than n columns. 



V. COMPLEXITY OF SIMULATING 
STABILIZER CIRCUITS 

The main benefit provided by the formalism of this pa- 
per above the existing techniques in the literature is that 
computing phases are effectively reduced to simple linear 
transformations, thereby simplifying the individual steps 
of simulating stabilizer circuits, which should reduce the 
burden of carrying out transformations in analytical in- 
vestigation and ad-hoc calculations. (This is not to say 
that these operations are asymptotically more efficient: 
we remark on this distinction in Section IVC1 ) However, 
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the reduction to linear algebra in itself also makes certain 
complexity theoretic results easier to prove, when consid- 
ering computational complexity classes which themselves 
are well characterized in linear algebraic terms. Using 
the techniques of Sections IIIII and IIV 1 we generalize the 
results of Aaronson and Gottesman [Xl| concerning the 
complexity of simulating stabilizer circuits to qudits of 
arbitrary dimension. Specifically, natural decision prob- 
lems concerning simulating stabilizer circuits on qudits 
of any fixed dimension d, involving at most a constant 
number of measurements, are complete for the complex- 
ity class coModrfL consisting of problems which are log- 
space reducible to determining whether a system of equa- 
tions mod d is feasible [30j . We describe these results in 
this section. 



A. Complexity of simulating unitary stabilizer 
circuits 

We define Definite Stabilizer Measurement to be 
the problem of deciding whether or not the state of a 
system of n qudits of dimension d, initially in a computa- 
tional state |q) and then acted on by a unitary stabilizer 
circuit, is stabilized by some particular P g Vf n . This is 
equivalent to the proposition that a measurement of the 
operator P would yield a record of "0" with certainty, 
being a +l-eigenstate of P. (Other possible outcomes 
h € "La may be considered instead by testing whether 
T -2hp stabilizes the state.) If the outcome "0" does not 
occur with certainty, it occurs either with probability 0, 
or probability 1/s (at most ^) for some integer s which 
divides d. 

In the special case of a P = —Z\ measurement on 
qubits (d = 2), Definite Stabilizer Measurement cor- 
responds to the problem Gottesman-Knill described by 
Aaronson and Gottesman 



11]. This problem belongs to 



the class ©L of problems which are log-space reducible 
to feasibility of systems of linear equations mod 2 and 
verifying coefficients of matrix products mod 2 (l2j ]. Re- 
lying on the result L® L = ©L [31(, Aaronson and Gottes- 
man describe an algorithm solving Gottesman-Knill on 
a logspace machine with access to an ©L oracle, which is 
used to repeatedly simulate initial segments of the stabi- 
lizer circuit in order to compute the effects of the phase 
corrections induced on intermediate states of the circuit. 
We show a more direct and generalized version of the 
result of Ref. [11] by characterizing the complexity of 
Definite Stabilizer Measurement. 

Just as ©L is the class of problems which are log- 
space reducible to verifying coefficients of matrix prod- 
ucts modulo 2, we may define the class coModdL as the 
class of decision problems which are log-space reducible 
to verifying coefficients of matrix products modulo d [l4| 
(see not e 1331 V In particular, this implies ©L = coMoc^l-. 
As Ref. [30] shows, testing whether a system of equations 
is feasible mod d is also a complete problem for coMod f /L 
for all d ^ 2; and in the particular case where d is a prime 



power, coMod^L may be characterized as those problems 
which are log-space reducible to evaluating (as opposed 
to verifying) coefficients of matrix products modulo d. 
These operations are in essence precisely what is required 
to simulate the transformations of a stabilizer tableau in 
a unitary stabilizer circuit, so that we may show: 

Theorem 7: For qudits of some fixed dimension d ^ 
2, Definite Stabilizer Measurement is complete for the 
class coModdL . 



Proof: As in the analysis of Section IIII D[ determining 
whether an n-qudit state, characterized by a (proper) 
stabilizer tableau 3}, is stabilized by an operator P = 
t~ 2S W p corresponds to determining whether there is a 
vector 5' © p' in the column-span of Tf such that P = 
T ~ 2S Wp>. This may be reduced to determining whether 
there is a solution to some system of equations At = 
(6 © p) (mod D) for some matrix A: 

• For d odd, in order for P = t~ 2S W p i, we require 
that p' - p = (mod d), and that 25' - 25 = 
]p', p] = (mod d), so that 6 © p is in the column 
span of Tf. We then set A :— Tf. 

• For d even, we still require p' — p = (mod d), 
but this is no longer sufficient to ensure p = p' € 
7jD\ and again we require 28' — 26 = [p',p] = 
(mod D). Then p' — p = dx for some x € {0, l} 2n , 
and 8' = 5 + | [p + dx, p] = | [x, p] (mod d). The 
vector c?x corresponds to some representation of the 
identity operator 1 = Wdx which accounts for the 
difference in the phase coefficients, by the formula 



-2<5 



' W = T -2S+[d*,v] W 



p+dx 



(47) 



by Lemma [5J As in Section IIII D[ we introduce 
auxiliary columns of the form Uo = d © and 
| [ej , p] © dej for 1 ^ j ^ 2n to the tableau 



u, 



Tf, chosen so that Uj + (S © p) also represents the 
operator P for each ^ j ^ 2n. Thus, it suf- 
fices to determine whether the system of equations 
At = [uo ui • • • U2n | Tf ] t = (6 © p) has solutions 
modulo 2c?. 



In each case, to test the feasibility of such a system of 
equations has solutions with a coModdL algorithm, it is 
not necessary to store Tf explicitly in the workspace; it 
suffices to be able to efficiently query individual coeffi- 
cients of Tf on demand, using only 0(log(n)) workspace. 

In the case where Tf represents a stabilizer state U |q) 
obtained by acting on a standard basis state |q) € %® n 
with a Clifford operator U oc Un ■ ■ ■ UiU\, we define Tf as 
the action of a sequence of operators Cn •■ - C%C\ acting 
on an initial tableau of the form T :— [q | ] , where 
each Cj is a conjugation tableau as in Eq. (|12p . We then 
consider two cases: 
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• Suppose d — p e for some prime p. Using the 
techniques of Ref. [30], for any prime power p e , 
we may compute the coefficients of such a ma- 
trix product modulo p e as part of a coModpL al- 
gorithm. Determining whether the system of equa- 
tions At = (6 © p) is feasible modulo each prime- 
power divisor p e is thus in coModpL. 

• For qudit dimensions having a prime power de- 
composition d — p^ ■■■p'j, 1 tor £ > 1, note that 
At = (<5fflp) (mod D) is feasible if and only if At = 
(5 © p) (mod p e ) is also feasible for every prime- 
power factor p e of D. Define the problem to 
be the problem Definite Stabilizer Measurement 
for a fixed qudit dimension d. Then we may char- 
acterize Lid as 

L d = L p? n L p? n • • • n L p7 . (48) 

By a standard normal form for coModpL [30L 
Prop. 3], it follows that e coModpL if and only 
if L p e g coModp»L (= coModpL) for each prime 
power p e which divides d. 

Thus, Definite Stabilizer Measurement for qudits of 
dimension d, is in coModpL for all integers d ^ 2. Finally, 
as simulating circuits consisting only of AX gates on qu- 
dits of dimension d is itself coMod^L-hard [33], it follows 
that Definite Stabilizer Measurement for dimension-d 
qudits is coMod^L-complete. □ 

We may generalize further, to consider the complexity 
of the problem of computing the output distribution of a 
measurement of a given observable. For a given qudit di- 
mension d and a probability distribution over Z^ , define 
Stabilizer Measurement to be the problem of deciding 
whether this distribution can be produced by measuring 
a given observable P = t~ 2S W p on the state produced 
by a specified unitary stabilizer circuit U acting on a 
given standard basis state |q). We restrict the distribu- 
tions taken as input to uniform distributions over a coset 
K + r/Zd C Zrf, where k, r\ € Z<j are parameters speci- 
fying the distribution. (The delta-peaked distributions 
of Definite Stabilizer Measurement correspond to the 
case 77 = d, and testing whether t~ 2k P stabilizes the 
state U |q); the uniform distribution over all of Z^ corre- 
sponds to 77 = 1, with k being redundant.) Consider 
the proper tableau Tf just prior to measurement: by 
Section fill D[ the outcome is uniformly distributed over 
K+r/Zd Q Z^ if and only if the row- vector <p '■= (0©p) T 2) 
has the property that rj = gcd(d, <f>i, 02, • • •)) an d the 
integer vector — S) © p] is generated modulo D 

by the columns of Tf (together with auxiliary columns 
Uj := I [e\j, p] © dej and Uo = d® in the case of d even, 
as in the proof of Theorem[7J). We may test both of these 
conditions by solving linear equations modulo D: this is 
obvious for the latter constraint, and we also have r\ = 
gcd(<i, 0i, 4>2, ■ ■ •) if and only if the system of equations 
i] = 4> T x (mod d) has solutions. As coModpL is closed 



under logical conjunctions, Stabilizer Measurement is 
thus also complete for coModpL. 

Note that Stabilizer Measurement is equivalent to 
determining whether the probability of obtaining a given 
outcome h is equal to 1/s for some s > 1, as this holds 
if and only if the outcome is uniformly distributed over 
h + (d/s)Zd- We therefore have: 

Theorem 8: For an initial state in the standard basis, 
any Pauli measurement observable P £ ■p®" j and any 
stabilizer circuit which performs a P measurement, the 
following problems are all coModrfL-complete: (a) verify- 
ing predictions of a deterministic measurement outcome 
h G Z<i, (b) verifying predictions of having probability p 
of obtaining an outcome h € Z^, and (c) verifying pre- 
dictions of the distribution of measurement outcomes. 

As a corollary, all of the problems described above may 
be simulated by 0(log(n) 2 )-depth boolean circuits, as 
coMod d L C NC 2 Q. 



B. The complexity of simulating stabilizer circuits 
with multiple measurements 

We now consider the way in which the results of the 
preceding section extend to complexity containments for 
stabilizer circuits with measurements, i.e. in which the 
evolution of the state under measurement must be ex- 
plicitly computed and where the outcomes may con- 
trol further operations. To simulate the transformation 
of a state under measurement — as opposed to deter- 
mining what the distribution of outcomes is, as in the 
Stabilizer Measurement problem — we must describe 
how a coModpL algorithm might carry out the calcula- 
tions described in Section HV CI 



1. Simulating evolution under a single measurement 

We first sketch an algorithm to verify any single co- 
efficient of an (extended) post-measurement tableau Tg' 
in coModpcL, for a qudit dimension of p e for some prime 
p, where the measurement outcome is somehow specified 
in the input and where we assume we may query coeffi- 
cients of the pre-measurement tableau T§. The tableau 
T§i is that which results from a Z r measurement acting 
on a pre-measurement tableau T§ via the procedure of 
Section IIV CI The motivation for restricting to prime- 
power qudit dimension is to describe a solution involving 
query access to T§ , which may be difficult to simulate in 
logarithmic space for composite d. We indicate how this 
extends to composite qudit dimensions d, in a manner 
similar to the proof of Theorem [71 before proceeding to 
the case of evolution under a sequence of measurements. 

To evaluate a coefficient of the post-measurement 
tableau T§' , we do not have to store the entire tableau as 
it is transformed, so long as we can efficiently reconstruct 



17 



the dependencies of the coefficient in question on the co- 
efficients of the pre-measurement tableau T§ on demand. 
We therefore describe how to reproduce which transfor- 
mations are performed on T§, and the impact of these 
transformations on the desired coefficient. 

The transformations performed on a tableau consist 
largely of column recombinations to clear the (n+r+l) st 
row of the tableau. These may be performed by invertible 
transformations, in which one column is added or sub- 
tracted from another some number of times; for a prime- 
power modulus, determining a combination in which the 
(n+r+l) st coefficient of one of the columns is sent to zero 
is easy, using comparisons and divisions of fixed-size in- 
tegers in {0,1,. — 1}. We may consider consecutive 
pairs of columns in turn — first considering combinations 
of the first column with the second column, then combi- 
nations of the second column with the third, and so on — 
and determine for each pair the transformation which will 
clear the (n+r+l) st coefficient of the left-most column in 
each case. By considering the effect of this sequence of 
column-combinations on the coefficients in other rows, we 
may determine how the coefficients in those rows trans- 
form, in order to determine what the value of any one 
given coefficient of T§> would be. We may do this us- 
ing enough workspace to store the (n+r+l) st coefficients 
of whichever two columns we consider at each step of 
the algorithm, as well as the coefficients for the same 
two columns in any other row we require; and enough 
workspace to carry out simple calculations, such as divi- 
sion, on fixed-width integers. Apart from the column re- 
combinations of the tableau, computing rj (and the scalar 
factor a by which we multiply the final non-zero column, 
as described following Eq. (02])) can be easily performed 
in constant space, as can s = d/rj for ij > 0. 

In the case of a power of an odd prime, the above 
suffices to determine all the coefficients of the tableau 
Tg which represents the same pre-measurement state as 
T§ , in which only a single generator Se fails to commute 
with Z r . In the case of d even, the phase correction 
block and the phase vector both involve row operations 
for each column operation. For a phase coefficient, we 
must also query off-diagonal coefficients from the phase 
correction block of the tableau, as it is being transformed; 
then some additional workspace is required to compute 
these coefficients. A phase correction block coefficient it- 
self may be subject to both column and row operations 
throughout the transformation of the stabilizer tableau: 
we may outline how these may be computed as follows. 
Suppose we wish to compute the value that the (h,j) 
coefficient of the phase correction block would have af- 
ter performing the fc th round of column recombinations, 
where 1 ^ h, j, k 2n. The phase correction block is an- 
tisymmetric by construction, throughout the transforma- 
tion of the tableau, which allows us to make the following 
observations: 

(i) If j = h, we return 0; and if j > h, we may instead 
compute the negation of the (j, h) coefficient of the 
phase correction block. 



(a) If j < h and k < h — 1, then none of the row- 
transformations on the phase correction block cor- 
responding to recombining the first k columns of 
the tableau have affected any coefficients in the h th 
row of the phase correction block, in which case we 
may simply simulate the effect of the first k column 
combinations on the (h,j) coefficient of the phase 
correction block, as for the Weyl block coefficients. 
(In particular, if k < j — 1, we may simply return 
the corresponding coefficient of T$.) 

(Hi) If j = k = h — 1, then the (h,j) coefficient is in 
principle affected by row-transformations between 
the (h — l) st row and the h th row. However, as 
the (h — 1, h— 1) coefficient of the phase correction 
block is zero, this row-operation has no effect, and 
we may reduce to the preceding case. 

(iv) If j < h ^ fc, the (h,j) coefficient of the phase 
block is affected by row-transformations which in 
general will have a non-trivial effect. We recur- 
sively compute coefficients b t j corresponding to the 
(t, j)-coefficient of the phase correction block after 
t column combinations, for j + 1 ^ t ^ h+1. Start- 
ing by computing foj+ij as in the preceding case, 
we compute each subsequent btj as follows: let b' t j 
be the value of the (t 1 , j)-coefficient of the phase 
correction block after the first t — 2 column com- 
binations, and then simulate the appropriate row- 
transformations with bt—i,j and b' t j to compute the 
value of btj after the (t — l) st column combination. 

All of the above can be performed in constant workspace, 
using at most two levels of recursive evaluation of the 
coefficients of the phase-correction block (as in the final 
case above). 

Having obtained a tableau Tg representing the pre- 
measurement group, for which only a single generator St 
fails to commute with Z r , we may easily describe the re- 
maining calculations required to determine coefficients of 
the tableau. The eigenvalue r 2t of the pre-measurement 
state with respect to the operator Z^ can be obtained 
by multiplying the measurement outcome specified for 
Z r at the input by s; we may then use a coMod p L ora- 
cle to solve Definite Stabilizer Measurement to deter- 
mine whether the specified outcome is possible, and if 
so, proceed with the computation (perhaps returning an 
error value otherwise). The column t © e r © £ describing 
the stabilizer arising from the measurement will be the 
(£ + l) st column of the tableau, for £ € {0, ^} as 
described in Eq. (j44j) . unless the column vector represent- 
ing the generator Sf can be expressed as a combination 
of t © e r © £ together with the other columns, modulo 
p e . We may determine this once more with a coMod p L 
oracle to solve systems of equations; and if there is indeed 
a solution, we omit the column corresponding to the old 
generator Si entirely. 

As a minor variation of the procedure of Section IIV CI 
we may switch the two columns representing the mea- 
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surement stabilizer r~ 2h Z r and the generator Sf, so that 
the former is always the £ th column and the latter the 
(£ + l) st column in the case that it is non-trivial. We may 
fixing the number of columns to the maximum of 2n as 
described in Section llV A 31 setting all additional columns 
to zero. All of the above suffice to compute any particular 
coefficient of the post-measurement tableau T§<, modulo 
p e ; if a test-value for the coefficient is provided as input, 
we may then test congruence modulo p e as well. 

Using the characterization of coModpL for arbitrary 
d ^ 2 (possibly divisible by multiple primes) described 
in Ref. [30, Prop. 3], this suffices to verify coefficients 
of a tableau transformed under measurement for arbi- 
trary qudit dimension as well. Any linear transforma- 
tions which are invertible modulo the prime-power di- 
visors pi 1 , j?2 2 , etc. of D are also invertible modulo D, 
by the Remainder Theorem; therefore each of the col- 
umn transformations modulo the prime powers p^ cor- 
respond to valid column transformations modulo D as 
well. By moving the column for the one generator of the 
group whose presence (or rather, whose status as a non- 
trivial generator) is uncertain to the (£ + l) st column, 
the columns in the tableaus for each prime-power divisor 
p e ? correspond to the same generators as one another; 
any prime-power divisor p e j for which some generator 
corresponds to the zero vector merely represents a Pauli 
stabilizer which is proportional to some Weyl operator 
W p ?j v . Thus, verifying the value of any given coefficient 
of a post-measurement tableau Tg' for qudits of any di- 
mension d 2 with respect to a Z r measurement is a 
problem contained in coModpL 



2. Simulating multiple measurements in arbitrary 
dimensions 

From the foregoing, it is straightforward to use known 
oracle-closure results for prime-power d [3l| to show that 
a stabilizer circuit involving any constant number of mea- 
surements may be simulated in coModpL, for qudits of 
arbitrary dimension d ^ 2. We proceed again along the 
same lines as the preceding section, by bounding the com- 
plexity for qudits of prime-power dimension, and then 
lifting to arbitrary dimension d. 

In the procedure above for simulating a measurement 
of a tableau in prime-power dimension p e , we assumed 
the ability to query individual coefficients of the tableau 
Tg which represents the pre-measurement state. Each 
coefficient of the post-measurement tableau may be eval- 
uated in coModpL, provided we supplement the computa- 
tion with an oracle for the coefficients of Tg . If Tg arises 
from the simulation of a stabilizer circuit — which may 
also be simulated in coMod p L — it follows that proper- 
ties of interest of the post-measurement tableau may be 
computed in coMod p L coM ° d f L . If Tg itself arises from a 
circuit which involves a single measurement, an oracle 
for evaluating coefficients of Tg can be implemented us- 
ing a coMod p L coM ° dl,L oracle, i.e. an oracle which itself 



has access to an oracle to evaluate the coefficients of the 
tableau immediately following the first measurement. 

Consider the problem of simulating a stabilizer circuit 
with k > 1 measurements, in the sense of computing co- 
efficients of the stabilizer tableau (possibly in order to de- 
termine whether a given outcome hk G occurs for the 
final measurement), given some sequence of intermediate 
measurement outcomes hi, ... , hk-i G Generalizing 
the description above, the simplest way to regard the 
computational complexity of this problem is to provide 
one layer of nested coMod p L oracles for each measure- 
ment, yielding a hierarchy of oracles for simulating each 
successive measurement. Simulating a stabilizer circuit 
with k measurements is then contained by the class 

coMod p L coMod P LCoModi>L 

with a tower of k— 1 oracles to simulate the first fc— 1 mea- 
surements. For any fixed k > 1 which is constant in the 
input size, this class is simply equal again to coMod p L, 
by the oracle closure results of Ref. [3l| for p prime (but 
see note [IH). Thus simulating the outcomes of any fixed 
number of measurements can be simulated in coMod p L 
for qudit dimension p e . 

For arbitrary dimensions, as we remark in the proof 
of Theorem the problem of simulating a stabilizer 
circuit on d dimensional qudits can be reduced to the 
corresponding problem for each prime-power factor p e 
of d. For instance, to determine whether a sequence of 
k G 0(1) measurements yields outcomes hi, . . . , hk G Z,j 
with certainty, we may simulate the same circuit with 
all coefficients evaluated mod p e , and also reduce each 
of the coefficients h t modulo p e . If instead we are inter- 
ested in whether a sequence of outcomes h t arise with 
some probability {e.g. so that hi occurs with probabil- 
ity pi, /12 occurs with probability p 2 |i given that hi oc- 
curred, etc.), we may express the probabilities in terms 
of a product of powers p n for each prime divisor p of 
d, and simulate for each measurement 1 ^ t ^ k whether 
or not the measurement outcome modulo p e occurs with 
probability p H by testing whether the measurement ob- 
servable t stabilizes the state at measurement. As 
these problems may be contained in coModpL for prime p, 
and characterize the yes instances for the same problems 
modulo d, these suffice to show: 

Theorem 9: For any constant k ^ 1, for any stabi- 
lizer state given as input (described either as a standard 
basis state or via an initial stabilizer tableau), for any 
sequence hi, . ■ ■ , hk G Z<j of measurement outcomes, and 
for any sequence of conditional probabilities pi, . . . ,pk of 
them occurring, the problem of determining whether the 
outcomes of the first k measurements are the outcomes 
hj for 1 ^ j ^ fc, each with probability P1P2 • • 'Pj, is 
coModrfL-complete. 

N.B. The issues described above relating to oracles do 
not arise if we are content to use a polynomial amount 
of workspace in the simulation. Simply using the tech- 
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niques of the preceding sections, we may easily simu- 
late arbitrary stabilizer circuits with any number of mea- 
surements on standard basis states in polynomial time, 
by transforming stabilizer tableaus which are explicitly 
stored in the work space. 



C. Run time complexity 

We conclude our considerations of the complexity of 
simulating stabilizer circuits with some remarks on more 
precise measures of complexity. 

The computational complexity bounds of Sections IV Al 
and EE do not impose any bounds whatsoever on the 
run-time complexity, except that it is polynomial by 
virtue of coMod d L C IMC 2 C P 14]. In practice, the 
degree of this polynomial will be quite large, even for the 
non-deterministic logspace Turing machines fl4| which 
are the usual model of computation used to define the 
class coModdL. Furthermore, using the techniques of 
Ref. [30| . the degree of the run-time required to solve 
systems of linear equations modulo d will increase with 
the value of d — and with the number of measurements 
being simulated — according to the size of the prime- 
power factors of d. However, this is an expected trade- 
off in time for the savings in the work-space required 
by the algorithms outlined in Sections IV Al and IV Bl As 
coModdL C NC 2 C DSPACE(log(n) 2 ), the large running- 
time may be regarded as a consequence of simulating sta- 
bilizer circuits with a constant number of rounds of mea- 
surement, but of any size, with either Turing machines 
which require only 0(log(n) 2 ) space or polynomial-size 
boolean circuits of depth only (9(log(n) 2 ). The tech- 
niques used to show containment in coModdL may be 
regarded as demonstrating upper bounds on the abstract 
computational power of stabilizer circuits (for any fixed 
number of measurement rounds) for any d ^ 2; a more 
practical approach to simulating stabilizer circuits, on 
a computer architecture having a small number of pro- 
cessors but a substantial amount of memory, is simply 
to maintain an explicit record of stabilizer tableaus and 
measurement outcomes. 

Apart from the results of Sections IV Al and IV Bl we 
may compare run-time complexity of the techniques of 
Sections IIIII and [TVl to those of Ref. [ill, Sec. Ill], for sim- 
ulating stabilizer circuits on qubits with the more tradi- 
tional "binary" representation on qubits (i.e. for the case 
d — 2 alone). In this setting, we consider the complexity 
only of simulating a generator of the Clifford group (such 
as a Pauli operator or a gate corresponding to one of 
the operators in Eq. (jT5j ) , or of single-qubit Z measure- 
ments. It must be noted that despite the elimination of 
quadratic phase corrections in the simulation of Clifford 
operations, there is no improvement in the asymptotic 
complexity of simulating a single Clifford group opera- 
tion compared to the results of Ref. [ll| ; it is 0(n) in 
each case for a single gate acting on a tableau represent- 
ing an n-qudit state (whether represented by a proper 



tableau, or an extended tableau). Furthermore, the al- 
gorithm in Ref. [Tl| for simulating measurements in the 
case d = 2 is more efficient than the algorithm presented 
in this article for simulating measurements, as the proce- 
dure presented in Section [VB I is essentially an extension 
of techniques of Ref. [l| to the case d ^ 2, for which 
Ref. [llT ] represents an improvement. Indeed, it seems 
likely that the techniques of Ref. [Ill Sec. Ill] could be 
easily extended using the linear formalism of this article 
to achieve a complexity of 0(n 2 ) for simulating Z mea- 
surements in any fixed prime dimension. However, as 
those techniques seem to rely on the fact that Z 2 (or Id 
for prime d) is a field, and in particular that the gener- 
ating set for any stabilizer group has size at most n, it 
is not immediately clear how such techniques would ex- 
tend even to the case of a prime power dimension. The 
improvement of the formalism of this article over that 
of Ref. [HI is not in terms of run-time complexity, but 
rather extending the space-bounded complexity theoretic 
results to arbitrary dimensions d ^ 2, and exploiting the 
reduction to linear algebra to do this more directly in the 
case of prime powers. 



VI. CONCLUSION 

We have presented techniques to simulate unitary sta- 
bilizer circuits on qudits of any constant dimension d ^ 2 
using linear transformations, and simulate terminal mea- 
surements by solving systems of linear equations, using 
Weyl operators to represent Pauli operators. In par- 
ticular: we demonstrate that the Clifford group can in 
each case be described effectively using only Pauli oper- 
ators and a group representation of the symplectic group 
Sp 2n (Z£)) over the integers modulo D, where D £ {d, 2d} 
is determined according to whether d is odd or even. We 
also present the first explicit treatment of the evolution of 
a state under measurement for composite dimension, and 
demonstrate how this may be achieved by linear transfor- 
mations of an extended tableau. This leads to a simple, 
easy to use formalism for simulating stabilizer circuits on 
qudits of arbitrary dimension. 

The motivation for the formalism of this article is not 
the actual run-time or circuit complexity, but rather a 
formalism for arbitrary dimension d ^ 2 which is as uni- 
form as possible, in which simulating individual opera- 
tions is as simple as possible (in the sense that each op- 
eration involves fewer natural arithmetic operations). In 
doing so, we reduce the problem of simulation substan- 
tially to standard techniques of linear algebra, which is il- 
lustrated by the directness of the proof of Theorem [7] that 
simulating the measurement of a unitary stabilizer circuit 
is complete for the class coModdL for any d ^ 2 (ex- 
tending the computational complexity results of Aaron- 
son and Gottesman [111]). 

There remain open questions with respect to the effi- 
cient simulation of stabilizer circuits, which tools of the 
sort presented in this article may help address. We have 
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shown that stabilizer circuits on qudits of any fixed di- 
mension, and any size, may be simulated in NC 2 provided 
it has a fixed number of measurements, as this task is 
complete for coModdL Is it possible to extend the num- 
ber of measurements, for instance, to 0(log(n)), or (if we 
restrict the controlled-Clifford operations to controlled- 
Paulis) to 0(n)l Can similar results be obtained if we 
keep the number of measurements fixed, and perhaps 
even fix the size of the circuit, but allow the qudit dimen- 
sion to be a prime power p e (of a fixed prime p) provided 
as input? Finally, what bounds can we obtain for simu- 
lating stabilizer circuits in a distributed classical compu- 
tational model based on linear transformations, such as 
linear network coding [35[? 
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Appendix A: Proper stabilizer tableaus and 
symplectic Clifford transformations 

We now show that any generating set for a stabilizer 
group over Hf n has a proper stabilizer tableau, for d 
prime or odd. We also show how the same techniques 
apply to prove the "lifting" claim of Theorem [5j 



1. Proper stabilizer tableaus for d prime and d odd 

For d either prime or odd, we wish to show that any 
sequence S\, . . . ,Sg G Vf n of commuting Pauli opera- 
tors (each with order at most d) may be represented 
by a sequence of vectors vi, . . . , G Zj^ and phases 
ipi, . . .(pi £ Zj3 such that Sj = t~ 2 ^ . , and for which 
[vh,vj] = (mod D) for all ^h,j ^ I. This is in fact 
trivially true in the case of d odd: from the hypothesis 
that the operators Sj commute pair-wise, it follows that 
[v/jjVj] = (mod D) by Lemma Q] for any Pauli vectors 
<fj®Vj representing the operators Sj as above. Thus, all 
stabilizer tableaus are proper for d odd. 

It remains to prove the result for the binary case d — 2. 
More generally, for any even d, we may construct a proper 
tableau for the operators Sj, provided S^S™ 2 ■ ■ ■ SJ 1 ' = 
1 only if each rrij is even. This holds, for example, for 
{Zi, Z2, • • • , Z n } on n qudits, or any family of operators 
Sj = UZjU^ obtained from them by conjugation by U G 
C n (d). We set Vi = vi, and for each subsequent j > 1, 
we construct Vj — Vj + dxj for Xj G Z|j" which satisfies 

[vj,Vj] = (mod 2d) and (Ala) 
[v h ,Vj] = (mod 2c?) for all 1 < h< j. (Alb) 



We may construct such a vector Vj as follows. The first 
constraint above ensures that = W Vj . As [v j 5 Vj] — 
0, we may simplify it to obtain [vj,dxj] = (mod 2d), 
or equivalently 

[vj,Xj-]=0 (mod 2). (A2a) 

For the second constraint above, note that as [v^, Vj] is a 
multiple of d for each 1 ^ h < j, we have [v/^Vj] = dbhj 
for some b^.j G Z; we may then expand this constraint 
to obtain 

= [v/i,Vj] = [vh,Vj] + [vh,dxj] 

= db hd + d[v h ,x j ] + 2d(i)[x h ,x j ] (mod 2d), (A2b) 

which we may further simplify to 

[v h ,Xj] = b h ,j (mod 2). (A2c) 

Let hj = bij §! + ■••+ &j and \Nj = \v% ■ ■ ■ Vj\- 

The constraints imposed by Eqs. (|A2ap and (|A2cp arc 
then equivalent to 

\Nja 2n = hj (mod 2) . (A3) 

Note that equations of the form S™ 1 • • • S™!^ SJ 1 = Al 
can only have solutions for A = 1 in stabilizer groups; and 
we have specifically ruled out the possibility that we can 
obtain such an expression for 1 using an odd exponent 
for Sj in our hypotheses. By Lemma[2l it follows that the 
vectors Vj are linearly independent modulo 2, in which 
case Eq. (| A3|) is a solvable system of linear equations over 
Z2 , with potentially many satisfactory solutions. Select- 
ing any one of them to fix a value of Vj , we then construct 
the next vector vy+i, and so forth until we have obtained 
a sequence of Weyl vectors vi, . . . , V£ forming a proper 
stabilizer tableau for the operators Si, . . . , St. 

As we note above, the above technique can be ap- 
plied for arbitrary qudit dimension d as well for stabilizer 
tableaus whose Weyl blocks have independent columns. 
For d composite, not all tableaus have this property. In 
particular, there exist stabilizer groups on n qudits which 
have more than n independent commuting generators, 
whose corresponding Weyl blocks therefore cannot have 
independent columns modulo 2 or any other prime. (See 
the beginning of Section IIV Al for an example, and Sec- 
tion IIV A 31 for a more general bound on the number of 
independent generators.) 

2. Symplectic transformations performed by 
Clifford operations 

We may apply similar techniques to the above to prove 
the "lifting" claim made in Theorem [SJ Note that if 
W Vj = UW^ J7 f for 1 < j < 2n, the operators W Vj will 
all have order d, and will be independent as well by virtue 
of the independence of the operators W$ . . In particular, 
the vectors Vj will be linearly independent modulo d. 
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We may then apply the same lifting technique as above, 
except that rather than the constraints of Eq. (|Aip . we 
impose the constraints 







(mod 2d) and 



(A4a) 



[vh,Vj] = S h<n+J ;- Sj.n+h (mod 2d) for h< j, (A4b) 

where 5 a ,b is the Kronecker delta. These same congru- 
ences will already hold modulo d (rather than 2d) by the 
preservation of commutation relations from the conjuga- 
tion by U. By constructing each Vj := Vj+dx.j as before, 
the same analysis may be applied to obtain constraints 
on each Xj sufficient to obtain the necessary inner prod- 
ucts modulo 2d as well. The matrix C = [vi • • • v 2n ] 
that we obtain as a result is symplectic modulo D. 



Appendix B: Supporting Lemmata concerning Weyl 
operators, stabilizer groups, and Pauli 
measurements 

We now present additional technical properties of Weyl 
operators and stabilizer groups which generalize the 
known results for d prime, to characterize of the evolu- 
tion of stabilizer states upon measurement. In particular, 
we show that this evolution is always uniformly random 
over a coset of the form K + rjZd for some k, r] G Z<j. Much 
of this appendix describes results and techniques which 
are both simple and standard for the case d = 2 (see 
e.g. Ref. [l|) and for d prime (Ref. Q). However, such 
techniques do not seem to have been explicitly described 
for d composite (for which technical obstacles exist, due 
e.g. to non-trivial operators Z a and X b which commute 
despite acting on a common qudit). 

We again adopt the convention described at the begin- 
ning of Section llTTDl that a measurement of a Pauli oper- 
ator P (of order at most d) stands for the measurement of 
some Hermitian operator H with spectral diameter less 
than 2ir, such that P = exp(iH), and that the outcome 
is some h G Id such that the result corresponds to the 
detection of a r 2a -eigenvector of P. We also suppose all 
measurements to be non-destructive, leaving a residual 
system which is in an eigenstate of the observable H, 
and thus of the unitary operator P. 



1. Spectral properties of Weyl operators 

The following simple results about Weyl operators are 
helpful to demonstrate how stabilizer groups transform 
under measurements in Sections IB 21 and IB 31 and may be 
of interest in the study of stabilizer codes over Z^. 

As W$ — 1 for all v G Z 2n , the eigenvalues of every 
Weyl operator are integer powers of r 2 . For d compos- 
ite, however, not all Weyl operators exhibit all possible 
powers of r 2 as eigenvalues. (This is true for d prime as 
well, if one includes the Weyl operator Wq = 1.) The 
eigenvalues that are exhibited by a Weyl operator W v , 



and their multiplicities, are governed in a simple way by 
the numerical relationships of the coefficients of v. 

Definition IX: For a vector v G Z^ 1 , the harmonic num- 
ber r](v) = gcd(i>i, . . . , v m , d) of v is the largest positive 
integer 77 ^ d such that v G T]7Jg. We will say that v is 
fundamental if r)(y) = 1. 

Accounting for the harmonic number is necessary to 
bridge the gap with the case of d prime 37] to obtain 
constructions for all d 2. Weyl operators W v with v 
fundamental have the greatest ability to "distinguish" be- 
tween different states; Weyl operators other than these 
have fewer eigenvalues, and each eigenvalue has higher 
multiplicity, so that their eigenspaces decompose as the 
sum of multiple eigenspaces of some other Weyl operator. 

Lemma 10: For each Weyl operator W v G Pf n , there 
exists a symplectic Clifford operator U G cC n (d) such 
that UW V W 



for any 1 ^ j ' ^ n. 



Proof: For v G Z 2n fundamental, using the decomposi- 
tion e.g. of Ref. [9|, Sec. IV], one may show that there 
is a symplectic operation C such that Cv = e n . This 
construction consists essentially of using row-reductions 
on a stabilizer tableau to compute the greatest common 
divisor of the coefficients of v; these row-reductions may 
be performed by left-multiplication by symplectic trans- 
formations. The operations needed to perform the row- 
reduction may be obtained by solving for coefficients aj 
such that a\Vi + ■ ■ ■ + ai n v?, n = 1. By the construction 
described in the proof of Lemma [6l there then exists a 
symplectic Clifford U G aC n {d) such that UW V W = Z n 
acting only on the n th qudit; we can then map this to 
any operator Zj by swaps. If v G Z 2n is not fundamen- 
tal, let u = v/r/(v); this vector is fundamental, so that 
there exists an operator U such that UW U W = Zj. Then 

uw v W = uw* v) W = zf v) . □ 

Note that the algorithm of Ref. [t| Sec. IV] yields the 
described unitary U in polynomial time, presented as a 
Clifford circuit of size 0(nlag(d)). 

Lemma 1 1 : For any n > and v G Z 2 ™ , the eigenvalues 
of Wv are all of the integer powers of e 27r "K v )/ d ; each of 
which occurs with multiplicity ?y(v) d n ~ l . 

Proof: By Lemma [TUl W v has the same spectrum as 
Zn . Taking Z n as an operator acting on Hf", its spec- 
trum consists of all of the integer powers of r 2 with mul- 
tiplicity d 71 ^ 1 (that is, equal multiplicity). As ?/(v) di- 
vides d, we may show that the spectrum of Z^ is all of 
the (d/rj(v)) th roots of unity also with equal multiplicity, 
which is to say ?y(v)cP -1 ; the same then holds for W v . □ 

Corollary 11a: For any v G Z 2 ^ 1 , the Weyl operator W v 
has order d/rj(y). In particular, it has order d if and only 
if v is fundamental. 
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2. Stabilizer groups of one-dimensional subspaces 
and commuting measurement observables 

In this section, we prove the usual connection of Pauli 
observables which commute with every element of the 
stabilizer group with that measurement yielding a de- 
terministic outcome, generalized to the setting of d 2 
arbitrary. To do so, we prove a characterization (which is 
well-known for d prime) of those stabilizer groups which 
describe unique pure states: they are maximal stabi- 
lizer groups, in the sense that any Pauli subgroup which 
strictly contains such a group is either nonabelian, or 
contains operators without +1 eigenvalues. (Proofs in 
the case of d prime typically make use of the fact that 
Z 2 ,™ is a field, which we cannot do for d composite.) 

Lemma 12: Let § = {Si, . . . , Sg} be a generating set 
for a stabilizer group Gg on n qudits. The following are 
equivalent: 

(i) Gg stabilizes a unique state; 

(ii) Gg is a maximal stabilizer group; 

(Hi) Gg is a stabilizer group of maximum size (and 
has cardinality d n ). 

We note property (Hi) above to emphasize the distinction 
from (ii), and for the sake of completeness. We will be 
interested primarily in certifying when a stabilizer group 
is inextensible, rather than when it has some particular 
cardinality. 

Proof: Let S C p® n be an arbitrary stabilizer group. 
For each generator Sj £ S, let 



d-i 



\Si\-l 



p=0 



(Bl) 



p=0 



where | Sj | is the multiplicative order of Sj . It is easy to 
show that this operator projects onto the +l-eigenspace 
of Sj. Then LTg = n^IIsj, • • • Hs e projects onto the joint 
+l-eigenspace of §. By expanding each of the projectors 
Hj , we may show that 



\G* 



(B2) 



seGs 



as the terms in the sum run over all distinct combina- 
tions of powers Sf 1 SJ? 2 ■ ■ ■ Sf l , generating each element 
of Gg. (As the group Gg is the direct product of the 
cyclic groups generated by the operators Sj, it follows 
that \Gg \ is equal to the product of the orders \Sj\.) 

It is easy to show that (i) <^=> (Hi): if S stabilizes a 
unique state, it follows that Tr(IIg) = 1, so that 



1 



Tr(n s ) 



1 



|G S 



E Tr (s) 



seG s 



Tr(l) 



(B3) 



where the final equality holds because 1 itself is the only 
element of Gg that has non-zero trace. Thus, \Gg\ = 



Tr(l) = d n ; and this cardinality is at a maximum, as 
Tr(II§) must be an integer. The converse is similar. 

It is also easy to show (i) =>• (ii). Suppose S stabi- 
lizes a unique state \ip), which is to say that Tr(II§) = 1. 
Let P £ Vf n be a Pauli operator which commutes with 
all of S. It is easy to show that § also stabilizes P 
by the uniqueness of l^), we therefore have P \ip) — A \tp) 
for some phase A. Then X _1 P stabilizes so that 



A^PILs 



A" 



|Gs| 



E PS 



seGs 



By hypothesis, we then have 
A- 1 



|Gg| 



E Tl i ps ) = Tr ( n s) = 1 



(B4) 



(B5) 
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As Tr (PS) ^ only if PS cx 1, this implies that there 
exists an operator S £ Gg such that S oc P^ 1 . As any 
two such operators would be proportional to one another, 
and can only be distinct if one of them failed to stabilize 
\ip), such an operator is unique. Similarly, as S oc AP _1 
both stabilize \tp), these are equal as well. Then either 
P £ Gg, in the case that A = 1; or the group obtained 
by extending Gg by P contains Al for A ^ 1, and is 
therefore not a stabilizer group. Thus the uniqueness of 
the state |^) entails that G§ is maximal as a stabilizer 
group in p® n . It remains to show that uniqueness of the 
stabilized state is a necessary condition for maximality 
as a stabilizer group. 

Suppose instead that Gg does not stabilize a unique 
state, and let |^o) ; IV'i) G i m g(IIs) be independent states 
stabilized by Gg. Consider the operator 



(B6) 



As the Weyl operators span the set of operators on H® n 
(by the corollary to Lemma on page ^ , there exists 
some operator W v such that 7 := Tr(W^r) 7^ 0. By 
construction, T has trace zero; and as TS = T = ST 
for every element S £ Gg, we have Tr(S^r) = for all 
S £ Gg as well. Thus W v is not proportional to any 
element of Gg. However, for any W s oc S £ Gg, we have 

7 = Tr^Wjr) = Tr^Wv^rW^)' 

= Tr(V 2 [ v < s JWvr) = t 2 [ v ' s 1 7 , (B7) 

where the penultimate equality holds by Lemma[T] Then 
[v, s] = (mod d) for all such s, so that S and W v com- 
mute. It follows that Wv commutes with all of Gg. 

While W v is not proportional to any element of Gg, 
there is a minimal integer 1 ^ s ^ d such that WJ is 
proportional to some P £ Gg . Let 77 = gcd(i>i, . . . V2n, d): 
by Lemma [Til W^ = W sv has eigenvalues consisting of 
integer powers of t 2sv with equal multiplicity. As P oc 
WJ has a non-trivial +l-eigenspace, the same is true for 
P. It follows that P = T 2sr,r W" for some ^ r < d. 
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Define S 1 * := r 2rir W v , which by construction satisfies 
Tr(Str) = r- 27 ' r Tr(Wjr) ^ 0, and each of whose inte- 
ger powers are either (a) an element of Gg or (b) not 
proportional to any element of Gg. In particular, as W v 
has order d/rj by construction, the only powers of S* 
which are proportional to the identity are in fact equal 
to 1. Consider then the group G obtained by extending 
Gg by 5"* . This group is abelian, Wv commutes 

with all of Gg • Because the intersection of (S*) and Gg 
is the subgroup (S%) by construction, we can decompose 
G into cosets of the form S^Gg for < t < s. Define the 
operator 



n 5 



G ^_ 

1 1 seG 



(B8) 



it is not difficult to show that IF 



rig , so that this is 



a projection. Furthermore, as S°l is the unique element 
of G proportional to the identity, we obtain 



Tr0) 
IGI 



>0, (B9) 



which implies that there exists a non-zero element of 
\ty) £ img(ng) which is a +l-eigenvector of each element 
of G. Thus G is a Pauli stabilizer group which strictly 
contains Gg, establishing (ii) =^> (i). □ 

This result has an important consequence for Pauli mea- 
surements: 

Corollary 12a: Let \ifj) £ Uf n be the unique +1- 
eigenstate of a stabilizer group Gg C Vf n , and P e Vf n 
an operator of order at most d which commutes with 
every operator in Gg. Then is undisturbed by mea- 
surements of P, and has a definite outcome h such that 

T -2hp g Gg _ 

Proof: By the preceding Lemma, there exists a scalar 
A such that XP 6 Gg . Then is an eigenvector of P 
and is undisturbed by measurement. As P has order at 
most d, we have P = t~ 2v W p for some integer vector 
p £ Z 2n : its eigenvalues consist of integer powers of of 
t 2 . The result of a P measurement on must be the 
h £ Z d for which T - 2h P e Gg. □ 



3. Measurement of observables which do not 
commute with all stabilizer generators 

Consider a maximal stabilizer group Gg C Vf n gen- 
erated by a set of operators § = {Si, . . . , Si}, stabilizing 
a state We consider the effect on \ip) of the mea- 
surement of an operator P oc W p such that P d = 1, but 
which does not commute with every Sj £ §. 



a. Reduction to the case of one generator not commuting 
with the measurement observable 



As in the case of d prime, we may reduce to the case 
where at most one stabilizer generator fails to commute 
with the measurement observable by considering different 
generating sets of Gg . For arbitrary d 2, we may do 
this as follows. 

For each 1 ^ j ^ £, consider coefficients ^ <fij < d 
(not all zero) such that PSjP^ = r 2<l,k S 3 . Note that 
the phases which are induced on Weyl operators W v by 
commutation with P depend on [p,v], which is linear 
in v: thus the phases induced by commutation with P on 
an arbitrary 5" € Gg by commutation with P is given by 
m\(f>i + • • • + mpjpz , where the integers rrij are exponents 
such that 5" = S™ 1 ■ ■ ■ S™ e . Let 77 = gcd(&, ...,&,<*): 
there exist integer vectors x e Z £+1 such that 



(BIO) 



From the above remarks, there then exists an element 
T = SI 1 --- S\ l £ Gg for which PTPt = t ^t. As 77 



divides each coefficient 



'.)■ 



let (fj — 



4>j/r]; we may then 
generate Gg by the operators S = {So, Si, . . . , Si}, set- 
ting S := T and Sj := SjT'W for 1 j < L 



The su. 



of stabilizers commuting with the 
measurement 



In the above construction, each generator Sj commutes 
with P for 1 ^ j ' ^ I. In the case that d is prime, 
they also generate the subgroup of Gg which commutes 
with the observable P. However, this does not hold for 
arbitrary d. For example, consider a Z\, measurement 
performed on a state \ip) ab of a system of two qudits a 
and b, where d — d\d2 (for d\,d% > 1) and where is 
stabilized by 



AX* F a |0)j0) 6 



(Bll) 



This state is stabilized by § = {Z- dl Z b , X a X dl }, as one 
may show by applying the transformations of Eq. p^|) 
to the operators {Z b , Z\} stabilizing the state |0) Q |0) b . 
Let Si = Z- d ^Z b and S 2 = X a X d \ Of the elements of §, 
only S% fails to commute with Z b ; and we obtain (f>\ = 0, 
4>2 = di . We may set z\ = and z 2 — 1 to obtain r\ = di , 



and from this define 



So 



: T 

Si 

5*2 



S°iSl 

S1T = Z- ai Z,,. 
S2T- 1 = 1 . 



x a x b 

— 



(B12) 



Note that {5i,if?2} alone fails to generate the operator 
X d 2 = s d 2 g q which commu tes with P. 
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In the general case, we may characterize the subgroup 
of G§ which commutes with P as follows. Consider an 
arbitrary operator S' = S™° S™ 1 ■ ■ ■ S™ e S G s which 
commutes with P. As 5™ J commutes with P for each 
1 ^ j ^ £, S™° must also commute with P, so that 

5™° = PS r na P^ = T 2r ' rn °S r na . (B13) 

Then mo is a multiple of dft], so that the operators 

§c = {Sq^ v , Si, §2, • ■ ■ , §i} generate the subgroup of Gs 
commuting with P. (In the familiar case of d prime, we 
necessarily have rj = 1, so that S^ 71 = 1 may be excluded 
as a generator; we then recover the result that S\, . . . , Si 
generate the commuting subgroup for d prime.) 



c. Characterization of the post-measurement state and 
distribution of outcomes 

Following the above analysis, suppose that we have 
a generating set § = {So, S%, . . . , S£\ in which only So 
fails to commute with P; specifically, we may suppose 
PSoP^ = t 2v So for some rj a divisor of d. From the pre- 
ceding section, §c = {S^ n ', Si, S2, ■ ■ • , Si} generates the 
subgroup of G§ which commutes with P. The operators 
of Sc also stabilize the post-measurement state: if we let 

d 

n h = -J2(T~ 2h py (Bi4) 

be the projector onto the r 2?1 -eigenspace of P, we have 

SIl h \i(>)=Il h S\ip)=Il h \iP) (B15) 

for any S £ Gg c . The post-measurement state is natu- 
rally also stabilized by R :— r~ 2h P = t~ 2u W p , where h 
is the actual measurement outcome. 

To consider the distribution of possible values of h G 
Zd, consider a Weyl operator W p oc P, and let u € Z f j 
be such that R = t~ 2u W p . By hypothesis, we have 
WpSoWl = t 2t >S , with W v S 3 Wl = Sj for all 1 < j < I 
Let s = djr\; just as Sq is the smallest power of So which 
commutes with P oc W p , we may show that 

S W;Sl = (r- 2d W p ) s = W£; (B16) 

it follows that W p commutes with all elements of §, and 
so is proportional to an element of G§ by Lemma [T2] 



Specifically, it is proportional to an element of Gg c , as 
Wp commutes with P. The post-measurement state is 
then stabilized both by t~ 2u W p and by W^. Thus the 
operator (t~ 2u W p ) s W p s = t~ 2su 1 stabilizes the post- 
measurement state, which implies that su = (mod d), 
i.e. that u is a multiple of d/s = r\. The result of mea- 
suring P then satisfies h e k + rfLd, where P = t~ 2k W p . 

The distribution of measurement outcomes within the 
coset k + r(Ld , and the relationships between them, 
straightforwardly generalize the standard results for d 
prime. Let Ilh again be the projection onto the r 2h - 
eigenspace of P as given in Eq. (|B14|) : then 

1 d 

s n h = -^Soir-^py 

d 

= - d Y.( T ~ 2h ~ 2Vp y S * = U h+r,S . (B17) 

It follows then that SoII^ = Ylh+ V \^)- By repeated 
application of So, we can iterate through all the ele- 
ments of k + rjZd , from which it follows that (a) all 
measurement outcomes in k + r\Ld are equally likely, as 
the vectors n K+??z have the same norm; and (b) the 
post-measurement state for the outcome h can be trans- 
formed to that for h + r\z by performing the operation 
Sq . Thus the "eliminated" stabilizer generator So of the 
pre-measurement state becomes a unitary byproduct op- 
erator (in the sense of Ref. 0]) which relates the possible 
post-measurement states; and the measurement outcome 
h is uniformly distributed over the coset k + rfE c i, where 
t 2k P is a Weyl operator. 

Note that the above analysis can accommodate the 
case of a deterministic measurement outcome, if we do 
not assume at the outset that P does not commute with 
all of Gs. If W p commutes with all elements of §, it fol- 
lows that rj = d, in which case the solution set to Eq. ([23)1 
is a singleton modulo d: the distribution of outcomes is 
then the "uniform distribution" on that singleton set, 
i.e. a delta distribution on the sole solution h. 

As we have remarked previously, the familiar case for 
d prime corresponds to rj = 1 provided that P does not 
commute with all elements of Gs. Then Eq. (|2"5|) triv- 
ializes, imposing no constraints on the value of h. The 
outcome is then uniformly distributed over all ft 6 Zj. 
The generalization for d composite is that the measure- 
ment outcome may in principle be uniformly distributed 
over a coset of any additive subgroup of Z^, and not just 
on some singleton set or on itself. 



